test.c 31 KB
Newer Older
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
4
 * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * 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 License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 * 
 ******************************************************************************/

Pedro Gonnet's avatar
Pedro Gonnet committed
21
22
/* Config parameters. */
#include "../config.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
23
24
25
26
27
28
29
30

/* Some standard headers. */
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <pthread.h>
#include <math.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
31
32
#include <float.h>
#include <limits.h>
33
#include <fenv.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
34
35
#include <omp.h>

Pedro Gonnet's avatar
Pedro Gonnet committed
36
/* Conditional headers. */
Pedro Gonnet's avatar
Pedro Gonnet committed
37
38
#ifdef HAVE_LIBZ
    #include <zlib.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
39
40
#endif

41
42
43
44
45
/* MPI headers. */
#ifdef WITH_MPI
    #include <mpi.h>
#endif

Pedro Gonnet's avatar
Pedro Gonnet committed
46
/* Local headers. */
47
#include "swift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
48
49
50
51
52
53

/* Ticks per second on this machine. */
#ifndef CPU_TPS
    #define CPU_TPS 2.67e9
#endif

54
55
56
57
58
/* Engine policy flags. */
#ifndef ENGINE_POLICY
    #define ENGINE_POLICY engine_policy_none
#endif

Pedro Gonnet's avatar
Pedro Gonnet committed
59
60
61
62
63
64
65

/**
 * @brief Mapping function to draw a specific cell (gnuplot).
 */

void map_cells_plot ( struct cell *c , void *data ) {

Pedro Gonnet's avatar
Pedro Gonnet committed
66
    int depth = *(int *)data;
Pedro Gonnet's avatar
Pedro Gonnet committed
67
68
    double *l = c->loc, *h = c->h;

Pedro Gonnet's avatar
Pedro Gonnet committed
69
    if ( c->depth <= depth ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
    
        printf( "%.16e %.16e %.16e\n" , l[0] , l[1] , l[2] );
        printf( "%.16e %.16e %.16e\n" , l[0]+h[0] , l[1] , l[2] );
        printf( "%.16e %.16e %.16e\n" , l[0]+h[0] , l[1]+h[1] , l[2] );
        printf( "%.16e %.16e %.16e\n\n\n" , l[0] , l[1]+h[1] , l[2] );
    
        printf( "%.16e %.16e %.16e\n" , l[0] , l[1] , l[2]+h[2] );
        printf( "%.16e %.16e %.16e\n" , l[0]+h[0] , l[1] , l[2]+h[2] );
        printf( "%.16e %.16e %.16e\n" , l[0]+h[0] , l[1]+h[1] , l[2]+h[2] );
        printf( "%.16e %.16e %.16e\n\n\n" , l[0] , l[1]+h[1] , l[2]+h[2] );
    
        printf( "%.16e %.16e %.16e\n" , l[0] , l[1] , l[2] );
        printf( "%.16e %.16e %.16e\n" , l[0] , l[1]+h[1] , l[2] );
        printf( "%.16e %.16e %.16e\n" , l[0] , l[1]+h[1] , l[2]+h[2] );
        printf( "%.16e %.16e %.16e\n\n\n" , l[0] , l[1] , l[2]+h[2] );
    
        printf( "%.16e %.16e %.16e\n" , l[0]+h[0] , l[1] , l[2] );
        printf( "%.16e %.16e %.16e\n" , l[0]+h[0] , l[1]+h[1] , l[2] );
        printf( "%.16e %.16e %.16e\n" , l[0]+h[0] , l[1]+h[1] , l[2]+h[2] );
        printf( "%.16e %.16e %.16e\n\n\n" , l[0]+h[0] , l[1] , l[2] +h[2]);
    
        printf( "%.16e %.16e %.16e\n" , l[0] , l[1] , l[2] );
        printf( "%.16e %.16e %.16e\n" , l[0] , l[1] , l[2]+h[2] );
        printf( "%.16e %.16e %.16e\n" , l[0]+h[0] , l[1] , l[2]+h[2] );
        printf( "%.16e %.16e %.16e\n\n\n" , l[0]+h[0] , l[1] , l[2] );
    
        printf( "%.16e %.16e %.16e\n" , l[0] , l[1]+h[1] , l[2] );
        printf( "%.16e %.16e %.16e\n" , l[0] , l[1]+h[1] , l[2]+h[2] );
        printf( "%.16e %.16e %.16e\n" , l[0]+h[0] , l[1]+h[1] , l[2]+h[2] );
        printf( "%.16e %.16e %.16e\n\n\n" , l[0]+h[0] , l[1]+h[1] , l[2] );
        
Pedro Gonnet's avatar
Pedro Gonnet committed
101
        if ( !c->split ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
102
            for ( int k = 0 ; k < c->count ; k++ )
Pedro Gonnet's avatar
Pedro Gonnet committed
103
104
105
106
                printf( "0 0 0 %.16e %.16e %.16e\n" ,
                    c->parts[k].x[0] , c->parts[k].x[1] , c->parts[k].x[2] );
            printf( "\n\n" );
            }
Pedro Gonnet's avatar
Pedro Gonnet committed
107
108
109
110
        /* else
            for ( int k = 0 ; k < 8 ; k++ )
                if ( c->progeny[k] != NULL )
                    map_cells_plot( c->progeny[k] , data ); */
Pedro Gonnet's avatar
Pedro Gonnet committed
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
    
        }

    }


/**
 * @brief Mapping function for checking if each part is in its box.
 */

/* void map_check ( struct part *p , struct cell *c , void *data ) {

    if ( p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] ||
         p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] ||
         p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] )
        printf( "map_check: particle %i is outside of its box.\n" , p->id );

    } */


131
132
133
134
/**
 * @brief Mapping function for neighbour count.
 */

Pedro Gonnet's avatar
Pedro Gonnet committed
135
136
137
138
139
void map_cellcheck ( struct cell *c , void *data ) {

    int k, *count = (int *)data;
    struct part *p;
    
140
141
    __sync_fetch_and_add( count , c->count );
    
Pedro Gonnet's avatar
Pedro Gonnet committed
142
143
144
145
146
    /* Loop over all parts and check if they are in the cell. */
    for ( k = 0 ; k < c->count ; k++ ) {
        p = &c->parts[k];
        if ( p->x[0] < c->loc[0] || p->x[1] < c->loc[1] || p->x[2] < c->loc[2] ||
             p->x[0] > c->loc[0] + c->h[0] || p->x[1] > c->loc[1] + c->h[1] || p->x[2] > c->loc[2] + c->h[2] ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
147
            printf( "map_cellcheck: particle at [ %.16e %.16e %.16e ] outside of cell [ %.16e %.16e %.16e ] - [ %.16e %.16e %.16e ].\n" ,
Pedro Gonnet's avatar
Pedro Gonnet committed
148
149
150
151
152
153
154
155
156
157
158
159
160
161
                p->x[0] , p->x[1] , p->x[2] ,
                c->loc[0] , c->loc[1] , c->loc[2] ,
                c->loc[0] + c->h[0] , c->loc[1] + c->h[1] , c->loc[2] + c->h[2] );
            error( "particle out of bounds!" );
            }
        }

    }


/**
 * @brief Mapping function for maxdepth cell count.
 */

162
163
164
165
166
167
168
169
170
171
172
173
174
void map_maxdepth ( struct cell *c , void *data ) {

    int maxdepth = ((int *)data)[0];
    int *count = &((int *)data)[1];
    
    // printf( "%e\n" , p->count );

    if ( c->depth == maxdepth )
        *count += 1;

    }


Pedro Gonnet's avatar
Pedro Gonnet committed
175
176
177
178
179
180
/**
 * @brief Mapping function for neighbour count.
 */

void map_count ( struct part *p , struct cell *c , void *data ) {

181
    double *wcount = (double *)data;
Pedro Gonnet's avatar
Pedro Gonnet committed
182
    
Pedro Gonnet's avatar
Pedro Gonnet committed
183
    // printf( "%i %e %e\n" , p->id , p->count , p->count_dh );
Pedro Gonnet's avatar
Pedro Gonnet committed
184

185
    *wcount += p->density.wcount;
Pedro Gonnet's avatar
Pedro Gonnet committed
186
187
188

    }

Pedro Gonnet's avatar
Pedro Gonnet committed
189
190
191
192
void map_wcount_min ( struct part *p , struct cell *c , void *data ) {

    struct part **p2 = (struct part **)data;
    
193
    if ( p->density.wcount < (*p2)->density.wcount )
Pedro Gonnet's avatar
Pedro Gonnet committed
194
195
196
197
198
199
200
201
        *p2 = p;

    }

void map_wcount_max ( struct part *p , struct cell *c , void *data ) {

    struct part **p2 = (struct part **)data;
    
202
    if ( p->density.wcount > (*p2)->density.wcount )
Pedro Gonnet's avatar
Pedro Gonnet committed
203
204
205
206
        *p2 = p;

    }

207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
void map_h_min ( struct part *p , struct cell *c , void *data ) {

    struct part **p2 = (struct part **)data;
    
    if ( p->h < (*p2)->h )
        *p2 = p;

    }

void map_h_max ( struct part *p , struct cell *c , void *data ) {

    struct part **p2 = (struct part **)data;
    
    if ( p->h > (*p2)->h )
        *p2 = p;

    }

Pedro Gonnet's avatar
Pedro Gonnet committed
225
226
227
228
229
230
231

/**
 * @brief Mapping function for neighbour count.
 */

void map_icount ( struct part *p , struct cell *c , void *data ) {

Pedro Gonnet's avatar
Pedro Gonnet committed
232
    // int *count = (int *)data;
Pedro Gonnet's avatar
Pedro Gonnet committed
233
234
235
    
    // printf( "%i\n" , p->icount );

Pedro Gonnet's avatar
Pedro Gonnet committed
236
    // *count += p->icount;
Pedro Gonnet's avatar
Pedro Gonnet committed
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263

    }


/**
 * @brief Mapping function to print the particle position.
 */

void map_dump ( struct part *p , struct cell *c , void *data ) {

    double *shift = (double *)data;

    printf( "%g\t%g\t%g\n" , p->x[0]-shift[0] , p->x[1]-shift[1] , p->x[2]-shift[2] );

    }


/**
 * @brief Read coordinates from a text file.
 *
 * @param fname The name of the coordinate file.
 * @param parts An array of #part in which to store the coordinates.
 * @param N The number of parts to read.
 */
 
void read_coords ( char *fname , struct part *parts , int N ) {

Pedro Gonnet's avatar
Pedro Gonnet committed
264
#ifdef HAVE_LIBZ
Pedro Gonnet's avatar
Pedro Gonnet committed
265
    gzFile fd;
Pedro Gonnet's avatar
Pedro Gonnet committed
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
    char buff[1024];
    int k;
    
    /* Open the given file. */
    if ( ( fd = gzopen( fname , "r" ) ) == NULL )
        error( "Failed to open coordinate file" );
        
    /* Read the coordinates into the part positions. */
    for ( k = 0 ; k < N ; k++ ) {
        if ( gzgets( fd , buff , 1024 ) == NULL )
            error( "Error reading coordinate file." );
        if ( sscanf( buff , "%lf %lf %lf" , &parts[k].x[0] , &parts[k].x[1] , &parts[k].x[2] ) != 3 ) {
            printf( "read_coords: failed to parse %ith entry.\n" , k );
            error( "Error parsing coordinate file." );
            }
        }
        
    /* Wrap it up. */
    gzclose( fd );
#else
Pedro Gonnet's avatar
Pedro Gonnet committed
286
287
288
289
290
291
292
293
294
295
296
297
298
299
    FILE *fd;
    int k;
    
    /* Open the given file. */
    if ( ( fd = fopen( fname , "r" ) ) == NULL )
        error( "Failed to open coordinate file" );
        
    /* Read the coordinates into the part positions. */
    for ( k = 0 ; k < N ; k++ ) {
        if ( fscanf( fd , "%lf %lf %lf" , &parts[k].x[0] , &parts[k].x[1] , &parts[k].x[2] ) != 3 ) {
            printf( "read_coords: failed to read %ith entry.\n" , k );
            error( "Error reading coordinate file." );
            }
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
300
301
302
303
        
    /* Wrap it up. */
    fclose( fd );
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
304
305
306
307
308
309
310
311
312
313
314
315
316
317

    }


/**
 * @brief Read cutoffs from a text file.
 *
 * @param fname The name of the cutoffs file.
 * @param parts An array of #part in which to store the cutoffs.
 * @param N The number of parts to read.
 */
 
void read_cutoffs ( char *fname , struct part *parts , int N ) {

Pedro Gonnet's avatar
Pedro Gonnet committed
318
#ifdef HAVE_LIBZ
Pedro Gonnet's avatar
Pedro Gonnet committed
319
    gzFile fd;
Pedro Gonnet's avatar
Pedro Gonnet committed
320
321
322
323
324
325
326
327
328
329
330
    char buff[1024];
    int k;
    
    /* Open the given file. */
    if ( ( fd = gzopen( fname , "r" ) ) == NULL )
        error( "Failed to open cutoff file" );
        
    /* Read the coordinates into the part positions. */
    for ( k = 0 ; k < N ; k++ ) {
        if ( gzgets( fd , buff , 1024 ) == NULL )
            error( "Error reading cutoff file." );
Pedro Gonnet's avatar
Pedro Gonnet committed
331
        if ( sscanf( buff , "%ef" , &parts[k].h ) != 1 ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
332
333
334
335
336
337
338
339
            printf( "read_cutoffs: failed to parse %ith entry.\n" , k );
            error( "Error parsing cutoff file." );
            }
        }
        
    /* Wrap it up. */
    gzclose( fd );
#else
Pedro Gonnet's avatar
Pedro Gonnet committed
340
341
342
343
344
345
346
347
348
    FILE *fd;
    int k;
    
    /* Open the given file. */
    if ( ( fd = fopen( fname , "r" ) ) == NULL )
        error( "Failed to open cutoff file" );
        
    /* Read the coordinates into the part positions. */
    for ( k = 0 ; k < N ; k++ ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
349
        if ( fscanf( fd , "%ef" , &parts[k].h ) != 1 ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
350
351
352
353
            printf( "read_cutoffs: failed to read %ith entry.\n" , k );
            error( "Error reading cutoff file." );
            }
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
354
355
356
357
        
    /* Wrap it up. */
    fclose( fd );
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
358
359
360
361
362
363
364
365
366
367
368
369
370
371

    }
    
    
/**
 * @brief Read id from a text file.
 *
 * @param fname The name of the id file.
 * @param parts An array of #part in which to store the dt.
 * @param N The number of parts to read.
 */
 
void read_id ( char *fname , struct part *parts , int N ) {

Pedro Gonnet's avatar
Pedro Gonnet committed
372
#ifdef HAVE_LIBZ
Pedro Gonnet's avatar
Pedro Gonnet committed
373
    gzFile fd;
Pedro Gonnet's avatar
Pedro Gonnet committed
374
375
376
377
378
379
380
381
382
383
384
    char buff[1024];
    int k;
    
    /* Open the given file. */
    if ( ( fd = gzopen( fname , "r" ) ) == NULL )
        error( "Failed to open id file" );
        
    /* Read the coordinates into the part positions. */
    for ( k = 0 ; k < N ; k++ ) {
        if ( gzgets( fd , buff , 1024 ) == NULL )
            error( "Error reading id file." );
385
        if ( sscanf( buff , "%lli" , &parts[k].id ) != 1 ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
386
387
388
389
390
391
392
393
            printf( "read_id: failed to parse %ith entry.\n" , k );
            error( "Error parsing id file." );
            }
        }
        
    /* Wrap it up. */
    gzclose( fd );
#else
Pedro Gonnet's avatar
Pedro Gonnet committed
394
395
396
397
398
399
400
401
402
    FILE *fd;
    int k;
    
    /* Open the given file. */
    if ( ( fd = fopen( fname , "r" ) ) == NULL )
        error( "Failed to open id file" );
        
    /* Read the coordinates into the part positions. */
    for ( k = 0 ; k < N ; k++ ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
403
        if ( fscanf( fd , "%lli" , &parts[k].id ) != 1 ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
404
405
406
407
408
            printf( "read_id: failed to read %ith entry.\n" , k );
            error( "Error reading id file." );
            }
        }

Pedro Gonnet's avatar
Pedro Gonnet committed
409
410
411
412
    /* Wrap it up. */
    fclose( fd );
#endif

Pedro Gonnet's avatar
Pedro Gonnet committed
413
414
415
416
417
418
419
420
421
422
423
424
425
    }
    
    
/**
 * @brief Read dt from a text file.
 *
 * @param fname The name of the dt file.
 * @param parts An array of #part in which to store the dt.
 * @param N The number of parts to read.
 */
 
void read_dt ( char *fname , struct part *parts , int N ) {

Pedro Gonnet's avatar
Pedro Gonnet committed
426
#ifdef HAVE_LIBZ
Pedro Gonnet's avatar
Pedro Gonnet committed
427
    gzFile fd;
Pedro Gonnet's avatar
Pedro Gonnet committed
428
429
430
431
432
433
434
435
436
437
438
    char buff[1024];
    int k;
    
    /* Open the given file. */
    if ( ( fd = gzopen( fname , "r" ) ) == NULL )
        error( "Failed to open dt file" );
        
    /* Read the coordinates into the part positions. */
    for ( k = 0 ; k < N ; k++ ) {
        if ( gzgets( fd , buff , 1024 ) == NULL )
            error( "Error reading id file." );
439
        if ( sscanf( buff , "%f" , &parts[k].dt ) != 1 )
Pedro Gonnet's avatar
Pedro Gonnet committed
440
441
442
443
444
445
            error( "Error parsing dt file." );
        }
        
    /* Wrap it up. */
    gzclose( fd );
#else
Pedro Gonnet's avatar
Pedro Gonnet committed
446
447
448
449
450
451
452
453
454
455
456
457
458
    FILE *fd;
    int k;
    
    /* Open the given file. */
    if ( ( fd = fopen( fname , "r" ) ) == NULL )
        error( "Failed to open dt file" );
        
    /* Read the coordinates into the part positions. */
    for ( k = 0 ; k < N ; k++ ) {
        if ( fscanf( fd , "%ef" , &parts[k].dt ) != 1 )
            error( "Error reading dt file." );
        }

Pedro Gonnet's avatar
Pedro Gonnet committed
459
460
461
462
    /* Wrap it up. */
    fclose( fd );
#endif

Pedro Gonnet's avatar
Pedro Gonnet committed
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
    }
    
    
/**
 * @brief Compute the average number of pairs per particle using
 *      a brute-force O(N^2) computation.
 *
 * @param dim The space dimensions.
 * @param parts The #part array.
 * @param N The number of parts.
 * @param periodic Periodic boundary conditions flag.
 */

void pairs_n2 ( double *dim , struct part *__restrict__ parts , int N , int periodic ) {

Pedro Gonnet's avatar
Pedro Gonnet committed
478
479
480
481
482
    int i, j, k, count = 0;
    // int mj, mk;
    // double maxratio = 1.0;
    double r2, dx[3], rho = 0.0;
    double rho_max = 0.0, rho_min = 100;
Pedro Gonnet's avatar
Pedro Gonnet committed
483
484
    
    /* Loop over all particle pairs. */
Pedro Gonnet's avatar
Pedro Gonnet committed
485
    #pragma omp parallel for schedule(dynamic), default(none), private(k,i,dx,r2), shared(periodic,parts,dim,N,stdout)
Pedro Gonnet's avatar
Pedro Gonnet committed
486
487
    for ( j = 0 ; j < N ; j++ ) {
        if ( j % 1000 == 0 ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
488
            printf( "pairs_n2: j=%i.\n" , j );
Pedro Gonnet's avatar
Pedro Gonnet committed
489
490
491
492
493
494
495
496
497
498
499
500
            fflush(stdout);
            }
        for ( k = j+1 ; k < N ; k++ ) {
            for ( i = 0 ; i < 3 ; i++ ) {
                dx[i] = parts[j].x[i] - parts[k].x[i];
                if ( periodic ) {
                    if ( dx[i] < -dim[i]/2 )
                        dx[i] += dim[i];
                    else if ( dx[i] > dim[i]/2 )
                        dx[i] -= dim[i];
                    }
                }
Pedro Gonnet's avatar
Pedro Gonnet committed
501
            r2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2];
Pedro Gonnet's avatar
Pedro Gonnet committed
502
            if ( r2 < parts[j].h*parts[j].h || r2 < parts[k].h*parts[k].h ) {
503
                runner_iact_density( r2 , NULL , parts[j].h , parts[k].h , &parts[j] , &parts[k] );
Pedro Gonnet's avatar
Pedro Gonnet committed
504
                /* if ( parts[j].h / parts[k].h > maxratio )
Pedro Gonnet's avatar
Pedro Gonnet committed
505
506
                    #pragma omp critical
                    {
Pedro Gonnet's avatar
Pedro Gonnet committed
507
                    maxratio = parts[j].h / parts[k].h;
Pedro Gonnet's avatar
Pedro Gonnet committed
508
509
                    mj = j; mk = k;
                    }
Pedro Gonnet's avatar
Pedro Gonnet committed
510
                else if ( parts[k].h / parts[j].h > maxratio )
Pedro Gonnet's avatar
Pedro Gonnet committed
511
512
                    #pragma omp critical
                    {
Pedro Gonnet's avatar
Pedro Gonnet committed
513
                    maxratio = parts[k].h / parts[j].h;
Pedro Gonnet's avatar
Pedro Gonnet committed
514
                    mj = j; mk = k;
Pedro Gonnet's avatar
Pedro Gonnet committed
515
                    } */
Pedro Gonnet's avatar
Pedro Gonnet committed
516
517
518
                }
            }
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
519
520
521
        
    /* Aggregate the results. */
    for ( k = 0 ; k < N ; k++ ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
522
        // count += parts[k].icount;
523
524
525
        rho += parts[k].density.wcount;
        rho_min = fmin( parts[k].density.wcount , rho_min );
        rho_min = fmax( parts[k].density.wcount , rho_max );
Pedro Gonnet's avatar
Pedro Gonnet committed
526
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
527
528
            
    /* Dump the result. */
Pedro Gonnet's avatar
Pedro Gonnet committed
529
    printf( "pairs_n2: avg. density per part is %.3f (nr. pairs %.3f).\n" , rho/N + 32.0/3 , ((double)count)/N );
Pedro Gonnet's avatar
Pedro Gonnet committed
530
531
    printf( "pairs_n2: densities are in [ %e , %e ].\n" , rho_min/N + 32.0/3 , rho_max/N + 32.0/3 );
    /* printf( "pairs_n2: maximum ratio between parts %i [%e,%e,%e] and %i [%e,%e,%e] is %.3f/%.3f\n" ,
Pedro Gonnet's avatar
Pedro Gonnet committed
532
533
        mj , parts[mj].x[0] , parts[mj].x[1] , parts[mj].x[2] ,
        mk , parts[mk].x[0] , parts[mk].x[1] , parts[mk].x[2] ,
Pedro Gonnet's avatar
Pedro Gonnet committed
534
535
536
537
538
539
        parts[mj].h , parts[mk].h ); fflush(stdout); */
    fflush(stdout);
            
    }


540
void pairs_single ( double *dim , long long int pid , struct part *__restrict__ parts , int N , int periodic ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
541
542
543
544
545

    int i, k;
    // int mj, mk;
    // double maxratio = 1.0;
    double r2, dx[3];
546
547
548
    float fdx[3];
    struct part p;
    // double ih = 12.0/6.25;
549
550
551
552
553
    
    /* Find "our" part. */
    for ( k = 0 ; k < N && parts[k].id != pid ; k++ );
    if ( k == N )
        error( "Part not found." );
554
    p = parts[k];
555
    printf( "pairs_single: part[%i].id == %lli.\n" , k , pid );
Pedro Gonnet's avatar
Pedro Gonnet committed
556
    
557
558
559
560
561
    p.rho = 0.0;
    p.density.wcount = 0.0;
    // p.icount = 0;
    p.rho_dh = 0.0;
            
Pedro Gonnet's avatar
Pedro Gonnet committed
562
563
    /* Loop over all particle pairs. */
    for ( k = 0 ; k < N ; k++ ) {
564
        if ( parts[k].id == p.id )
Pedro Gonnet's avatar
Pedro Gonnet committed
565
566
            continue;
        for ( i = 0 ; i < 3 ; i++ ) {
567
            dx[i] = p.x[i] - parts[k].x[i];
Pedro Gonnet's avatar
Pedro Gonnet committed
568
569
570
571
572
573
            if ( periodic ) {
                if ( dx[i] < -dim[i]/2 )
                    dx[i] += dim[i];
                else if ( dx[i] > dim[i]/2 )
                    dx[i] -= dim[i];
                }
574
            fdx[i] = dx[i];
Pedro Gonnet's avatar
Pedro Gonnet committed
575
            }
576
577
578
579
580
        r2 = fdx[0]*fdx[0] + fdx[1]*fdx[1] + fdx[2]*fdx[2];
        if ( r2 < p.h*p.h ) {
            runner_iact_nonsym_density( r2 , fdx , p.h , parts[k].h , &p , &parts[k] );
            /* printf( "pairs_simple: interacting particles %lli [%i,%i,%i] and %lli [%i,%i,%i], r=%e.\n" ,
                pid , (int)(p.x[0]*ih) , (int)(p.x[1]*ih) , (int)(p.x[2]*ih) ,
581
                parts[k].id , (int)(parts[k].x[0]*ih) , (int)(parts[k].x[1]*ih) , (int)(parts[k].x[2]*ih) ,
582
                sqrtf(r2) ); */
Pedro Gonnet's avatar
Pedro Gonnet committed
583
584
585
586
            }
        }
        
    /* Dump the result. */
587
    printf( "pairs_single: wcount of part %lli (h=%e) is %f.\n" , p.id , p.h , p.density.wcount + 32.0/3 );
Pedro Gonnet's avatar
Pedro Gonnet committed
588
    fflush(stdout);
Pedro Gonnet's avatar
Pedro Gonnet committed
589
    
Pedro Gonnet's avatar
Pedro Gonnet committed
590
591
592
593
594
595
596
597
598
599
600
601
602
    }


/**
 * @brief Find the pairs of a single particle
 *
 * @param dim The space dimensions.
 * @param parts The #part array.
 * @param N The number of parts.
 * @param periodic Periodic boundary conditions flag.
 * @param target the index of the target particle.
 */

Pedro Gonnet's avatar
Pedro Gonnet committed
603
void pairs_single_old ( double *dim , struct part *__restrict__ parts , int N , int periodic , int target ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
604
605

    int i, k, tid;
Pedro Gonnet's avatar
Pedro Gonnet committed
606
    double r, tx[3], th, dx[3];
Pedro Gonnet's avatar
Pedro Gonnet committed
607
608
609
610
    
    /* Get the target position and radius. */
    for ( k = 0 ; k < 3 ; k++ )
        tx[k] = parts[target].x[k];
Pedro Gonnet's avatar
Pedro Gonnet committed
611
    th = parts[target].h;
Pedro Gonnet's avatar
Pedro Gonnet committed
612
613
614
    tid = parts[target].id;
    
    /* Loop over all particle pairs. */
Pedro Gonnet's avatar
Pedro Gonnet committed
615
    #pragma omp parallel for schedule(dynamic), default(none), private(k,i,dx,r), shared(target,tx,th,tid,periodic,parts,dim,N)
Pedro Gonnet's avatar
Pedro Gonnet committed
616
617
618
619
620
621
622
623
624
625
626
627
628
    for ( k = 0 ; k < N ; k++ ) {
        if ( k == target )
            continue;
        for ( i = 0 ; i < 3 ; i++ ) {
            dx[i] = tx[i] - parts[k].x[i];
            if ( periodic ) {
                if ( dx[i] < -dim[i]/2 )
                    dx[i] += dim[i];
                else if ( dx[i] > dim[i]/2 )
                    dx[i] -= dim[i];
                }
            }
        r = sqrt( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
Pedro Gonnet's avatar
Pedro Gonnet committed
629
        if ( r < th )
630
            printf( "pairs_single: %i %lli [%e,%e,%e] %e\n" ,
Pedro Gonnet's avatar
Pedro Gonnet committed
631
632
633
634
                tid , parts[k].id , dx[0] , dx[1] , dx[2] , r );
        }
            
    }
635
636
637
638
639
640
641
642
643
644
645
646
    
    
/**
 * @brief Test the kernel function by dumping it in the interval [0,1].
 *
 * @param N number of intervals in [0,1].
 */
 
void kernel_dump ( int N ) {

    int k;
    float x, w, dw_dx;
647
648
649
    float x4[4] = {0.0f,0.0f,0.0f,0.0f};
    float w4[4] = {0.0f,0.0f,0.0f,0.0f};
    // float dw_dx4[4] __attribute__ ((aligned (16)));
650
651

    for ( k = 0 ; k <= N ; k++ ) {
652
        x = ((float)k) / N;
Pedro Gonnet's avatar
Pedro Gonnet committed
653
        x4[3] = x4[2]; x4[2] = x4[1]; x4[1] = x4[0]; x4[0] = x;
654
        kernel_deval( x , &w , &dw_dx );
655
        // kernel_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
Pedro Gonnet's avatar
Pedro Gonnet committed
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
        printf( " %e %e %e %e %e %e %e\n" , x , w , dw_dx , w4[0] , w4[1] , w4[2] , w4[3] );
        }

    }


/**
 * @brief Test the density function by dumping it for two random parts.
 *
 * @param N number of intervals in [0,1].
 */
 
void density_dump ( int N ) {

    int k;
    float r2[4] = {0.0f,0.0f,0.0f,0.0f}, hi[4], hj[4];
    struct part *pi[4], *pj[4], Pi[4], Pj[4];
    
    /* Init the interaction parameters. */
    for ( k = 0 ; k < 4 ; k++ ) {
676
677
        Pi[k].mass = 1.0f; Pi[k].rho = 0.0f; Pi[k].density.wcount = 0.0f;
        Pj[k].mass = 1.0f; Pj[k].rho = 0.0f; Pj[k].density.wcount = 0.0f;
Pedro Gonnet's avatar
Pedro Gonnet committed
678
679
680
681
682
683
684
685
686
        hi[k] = 1.0;
        hj[k] = 1.0;
        pi[k] = &Pi[k];
        pj[k] = &Pj[k];
        }

    for ( k = 0 ; k <= N ; k++ ) {
        r2[3] = r2[2]; r2[2] = r2[1]; r2[1] = r2[0];
        r2[0] = ((float)k) / N;
687
        Pi[0].density.wcount = 0; Pj[0].density.wcount = 0;
Pedro Gonnet's avatar
Pedro Gonnet committed
688
        runner_iact_density( r2[0] , NULL , hi[0] , hj[0] , &Pi[0] , &Pj[0] );
689
690
691
692
693
        printf( " %e %e %e" , r2[0] , Pi[0].density.wcount , Pj[0].density.wcount );
        Pi[0].density.wcount = 0; Pj[0].density.wcount = 0;
        Pi[1].density.wcount = 0; Pj[1].density.wcount = 0;
        Pi[2].density.wcount = 0; Pj[2].density.wcount = 0;
        Pi[3].density.wcount = 0; Pj[3].density.wcount = 0;
Pedro Gonnet's avatar
Pedro Gonnet committed
694
        runner_iact_vec_density( r2 , NULL , hi , hj , pi , pj );
695
        printf( " %e %e %e %e\n" , Pi[0].density.wcount , Pi[1].density.wcount , Pi[2].density.wcount , Pi[3].density.wcount );
696
697
698
        }

    }
Pedro Gonnet's avatar
Pedro Gonnet committed
699
700
701
702
703
704
705
706
707


/**
 * @brief Main routine that loads a few particles and generates some output.
 *
 */
 
int main ( int argc , char *argv[] ) {

708
    int c, icount, j, k, N = 100, periodic = 1;
Pedro Gonnet's avatar
Pedro Gonnet committed
709
    int nr_threads = 1, nr_queues = -1, runs = INT_MAX;
710
    int data[2];
Pedro Gonnet's avatar
Pedro Gonnet committed
711
    double dim[3] = { 1.0 , 1.0 , 1.0 }, shift[3] = { 0.0 , 0.0 , 0.0 };
712
    double h_max = -1.0 , scaling = 1.0;
Pedro Gonnet's avatar
Pedro Gonnet committed
713
    double clock = DBL_MAX;
714
    struct part *parts = NULL;
Pedro Gonnet's avatar
Pedro Gonnet committed
715
    struct space s;
Pedro Gonnet's avatar
Pedro Gonnet committed
716
    struct engine e;
717
    char ICfileName[200];
718
    float dt_max = 0.0f;
Pedro Gonnet's avatar
Pedro Gonnet committed
719
    ticks tic;
720
    int nr_nodes = 1, myrank = 0, grid[3] = { 1 , 1 , 1 };
Pedro Gonnet's avatar
Pedro Gonnet committed
721
    
722
723
724
    /* Choke on FP-exceptions. */
    feenableexcept( FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW );
    
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
#ifdef WITH_MPI
    /* Start by initializing MPI. */
    int res, prov;
    if ( ( res = MPI_Init_thread( &argc , &argv , MPI_THREAD_MULTIPLE , &prov ) ) != MPI_SUCCESS )
        error( "Call to MPI_Init failed with error %i." , res );
    if ( prov != MPI_THREAD_MULTIPLE )
        error( "MPI does not provide the level of threading required (MPI_THREAD_MULTIPLE)." );
    if ( ( res = MPI_Comm_size( MPI_COMM_WORLD , &nr_nodes ) != MPI_SUCCESS ) )
        error( "MPI_Comm_size failed with error %i." , res );
    if ( ( res = MPI_Comm_rank( MPI_COMM_WORLD , &myrank ) ) != MPI_SUCCESS )
        error( "Call to MPI_Comm_rank failed with error %i." , res );
    if ( ( res = MPI_Comm_set_errhandler( MPI_COMM_WORLD , MPI_ERRORS_RETURN ) ) != MPI_SUCCESS )
        error( "Call to MPI_Comm_set_errhandler failed with error %i." , res );
    if ( myrank == 0 )
        message( "MPI is up and running with %i nodes." , nr_nodes );
    fflush(stdout);
#endif
    
    
Pedro Gonnet's avatar
Pedro Gonnet committed
744
745
    /* Init the space. */
    bzero( &s , sizeof(struct space) );
746
747

    /* Parse the options */
748
    while ( ( c = getopt( argc , argv  , "a:c:d:f:g:m:q:r:s:t:w:z:" ) ) != -1 )
749
750
751
752
753
      switch( c )
	{
	case 'a':
	  if ( sscanf( optarg , "%lf" , &scaling ) != 1 )
	    error( "Error parsing cutoff scaling." );
754
755
	  if ( myrank == 0 )
        message( "scaling cutoff by %.3f." , scaling ); fflush(stdout);
756
	  break;
Pedro Gonnet's avatar
Pedro Gonnet committed
757
758
759
	case 'c':
	  if ( sscanf( optarg , "%lf" , &clock ) != 1 )
	    error( "Error parsing clock." );
760
761
	  if ( myrank == 0 )
        message( "clock set to %.3e." , clock ); fflush(stdout);
Pedro Gonnet's avatar
Pedro Gonnet committed
762
	  break;
763
	case 'd':
764
	  if ( sscanf( optarg , "%f" , &dt_max ) != 1 )
765
	    error( "Error parsing timestep." );
766
767
	  if ( myrank == 0 )
        message( "dt set to %e." , dt_max ); fflush(stdout);
768
769
770
771
	  break;
	case 'f':
	  if( !strcpy(ICfileName, optarg))
	    error("Error parsing IC file name.");
772
773
774
775
776
777
778
779
	  if ( myrank == 0 )
        message("IC to be read from file '%s'.", ICfileName);
	  break;
	case 'g':
	  if ( sscanf( optarg , "%i %i %i" , &grid[0] , &grid[1] , &grid[2] ) != 3 )
	    error( "Error parsing grid." );
	  if ( myrank == 0 )
        message( "grid set to [ %i %i %i ]." , grid[0] , grid[1] , grid[2] ); fflush(stdout);
780
	  break;
781
782
783
	case 'm':
	  if ( sscanf( optarg , "%lf" , &h_max ) != 1 )
	    error( "Error parsing h_max." );
784
785
	  if ( myrank == 0 )
        message( "maximum h set to %e." , h_max ); fflush(stdout);
786
787
788
789
790
791
792
793
794
795
796
797
	  break;
	case 'q':
	  if ( sscanf( optarg , "%d" , &nr_queues ) != 1 )
	    error( "Error parsing number of queues." );
	  break;
	case 'r':
	  if ( sscanf( optarg , "%d" , &runs ) != 1 )
	    error( "Error parsing number of runs." );
	  break;
	case 's':
	  if ( sscanf( optarg , "%lf %lf %lf" , &shift[0] , &shift[1] , &shift[2] ) != 3 )
	    error( "Error parsing shift." );
798
799
	  if ( myrank == 0 )
        message( "will shift parts by [ %.3f %.3f %.3f ]." , shift[0] , shift[1] , shift[2] );
800
801
802
803
804
805
	  break;
	case 't':
	  if ( sscanf( optarg , "%d" , &nr_threads ) != 1 )
	    error( "Error parsing number of threads." );
	  omp_set_num_threads( nr_threads );
	  break;
806
807
808
	case 'w':
	  if ( sscanf( optarg , "%d" , &space_subsize ) != 1 )
	    error( "Error parsing sub size." );
809
810
	  if ( myrank == 0 )
        message( "sub size set to %i." , space_subsize );
811
	  break;
812
813
814
	case 'z':
	  if ( sscanf( optarg , "%d" , &space_splitsize ) != 1 )
	    error( "Error parsing split size." );
815
816
	  if ( myrank == 0 )
        message( "split size set to %i." , space_splitsize );
817
818
819
820
821
822
	  break;
	case '?':
	  error( "Unknown option." );
	  break;

	}
823
            
824

825
    /* How large are the parts? */
826
827
    if ( myrank == 0 )
        message( "sizeof(struct part) is %li bytes." , (long int)sizeof( struct part ));
828
829
830

    /* Read particles and space information from (GADGET) IC */
    tic = getticks();
831
    read_ic( ICfileName , dim , &parts , &N , &periodic );
832
833
    if ( myrank == 0 )
        message( "reading particle properties took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
834
835
    
    /* Apply h scaling */
Pedro Gonnet's avatar
Pedro Gonnet committed
836
    if( scaling != 1.0 )
837
      for ( k = 0 ; k < N ; k++ )
838
	    parts[k].h *= scaling;
839
840
841
842
    
    /* Apply shift */
    if(shift[0] !=0 || shift[1] !=0 || shift[2] !=0 )
      for ( k = 0 ; k < N ; k++ ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
843
844
845
	    parts[k].x[0] += shift[0];
	    parts[k].x[1] += shift[1];
	    parts[k].x[2] += shift[2];
846
847
      }

Pedro Gonnet's avatar
Pedro Gonnet committed
848
849
850
851
852
853
854
    /* Set default number of queues. */
    if ( nr_queues < 0 )
        nr_queues = nr_threads;
            
    /* Initialize the space with this data. */
    tic = getticks();
    space_init( &s , dim , parts , N , periodic , h_max );
855
856
    if ( myrank == 0 )
        message( "space_init took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
857

858
    /* Set the default time step to 1.0f. */
859
860
    if ( myrank == 0 )
        message( "dt_max is %e." , dt_max );
861
    
Pedro Gonnet's avatar
Pedro Gonnet committed
862
    /* Say a few nice things about the space we just created. */
863
864
865
866
867
868
869
870
    if ( myrank == 0 ) {
        message( "space dimensions are [ %.3f %.3f %.3f ]." , s.dim[0] , s.dim[1] , s.dim[2] );
        message( "space %s periodic." , s.periodic ? "is" : "isn't" );
        message( "highest-level cell dimensions are [ %i %i %i ]." , s.cdim[0] , s.cdim[1] , s.cdim[2] );
        message( "%i parts in %i cells." , s.nr_parts , s.tot_cells );
        message( "maximum depth is %d." , s.maxdepth );
        // message( "cutoffs in [ %g %g ]." , s.h_min , s.h_max ); fflush(stdout);
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
871
    
Pedro Gonnet's avatar
Pedro Gonnet committed
872
    /* Verify that each particle is in it's propper cell. */
873
874
875
876
877
    if ( myrank == 0 ) {
        icount = 0;
        space_map_cells_pre( &s , 0 , &map_cellcheck , &icount );
        message( "map_cellcheck picked up %i parts." , icount );
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
878
    
879
880
881
882
883
    if ( myrank == 0 ) {
        data[0] = s.maxdepth; data[1] = 0;
        space_map_cells_pre( &s , 0 , &map_maxdepth , data );
        message( "nr of cells at depth %i is %i." , data[0] , data[1] );
        }
884
    
Pedro Gonnet's avatar
Pedro Gonnet committed
885
886
887
    /* Dump the particle positions. */
    // space_map_parts( &s , &map_dump , shift );
    
888
889
    
    /* Initialize the engine with this space. */
Pedro Gonnet's avatar
Pedro Gonnet committed
890
    tic = getticks();
891
892
893
894
895
896
897
898
899
900
    engine_init( &e , &s , dt_max , nr_threads , nr_queues , nr_nodes , myrank , ENGINE_POLICY | engine_policy_steal );
    if ( myrank == 0 )
        message( "engine_init took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);

#ifdef WITH_MPI
    /* Split the space. */
    if ( nr_nodes != grid[0]*grid[1]*grid[2] )
        error( "Grid size does not match number of nodes." );
    engine_split( &e , grid );
#endif
901
902

    /* Write the state of the system as it is before starting time integration. */
903
904
905
906
907
    if ( myrank == 0 ) {
        tic = getticks();
        write_output(&e);
        message( "writing particle properties took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
908
    
909
910
911
912
913
914
    /* Init the runner history. */
    #ifdef HIST
    for ( k = 0 ; k < runner_hist_N ; k++ )
        runner_hist_bins[k] = 0;
    #endif
    
Pedro Gonnet's avatar
Pedro Gonnet committed
915
916
    /* Inauguration speech. */
    if ( runs < INT_MAX )
917
        message( "starting for %i steps with %i threads and %i queues..." , runs , e.nr_threads , e.sched.nr_queues );
Pedro Gonnet's avatar
Pedro Gonnet committed
918
    else
919
        message( "starting for t=%.3e with %i threads and %i queues..." , clock , e.nr_threads , e.sched.nr_queues );
Pedro Gonnet's avatar
Pedro Gonnet committed
920
921
922
    fflush(stdout);
    
    /* Legend. */
923
924
    if ( myrank == 0 )
        printf( "# step time e_tot e_kin e_temp dt dt_step count dt_min dt_max\n" );
Pedro Gonnet's avatar
Pedro Gonnet committed
925
    
Pedro Gonnet's avatar
Pedro Gonnet committed
926
    /* Let loose a runner on the space. */
Pedro Gonnet's avatar
Pedro Gonnet committed
927
    for ( j = 0 ; j < runs && e.time < clock ; j++ ) {
928
    
929
        // message( "starting run %i/%i (t=%.3e) with %i threads and %i queues..." , j+1 , runs , e.time , e.nr_threads , e.nr_queues ); fflush(stdout);
930
        timers_reset( timers_mask_all );
Pedro Gonnet's avatar
Pedro Gonnet committed
931
932
933
934
        #ifdef COUNTER
            for ( k = 0 ; k < runner_counter_count ; k++ )
                runner_counter[k] = 0;
        #endif
935
936
        
        /* Take a step. */
Pedro Gonnet's avatar
Pedro Gonnet committed
937
        engine_step( &e );
Pedro Gonnet's avatar
Pedro Gonnet committed
938
        
939
940
941
        if ( myrank == 0 && j % 100 == 0 )
            write_output(&e);
                
Pedro Gonnet's avatar
Pedro Gonnet committed
942
        /* Dump a line of agregate output. */
943
944
945
946
947
948
949
950
951
952
        if ( myrank == 0 ) {
            printf( "%i %e %.16e %.16e %.16e %.3e %.3e %i %.3e %.3e" ,
                j , e.time ,
                e.ekin+e.epot , e.ekin , e.epot ,
                e.dt , e.dt_step , e.count_step ,
                e.dt_min , e.dt_max );
            for ( k = 0 ; k < timer_count ; k++ )
                printf( " %.3f" , ((double)timers[k])/CPU_TPS*1000 );
            printf( "\n" ); fflush(stdout);
            }
953
        
Pedro Gonnet's avatar
Pedro Gonnet committed
954
955
        }
        
956
957
958
959
960
961
962
963
964
965
    /* Print the values of the runner histogram. */
    #ifdef HIST
        printf( "main: runner histogram data:\n" );
        for ( k = 0 ; k < runner_hist_N ; k++ )
            printf( " %e %e %e\n" ,
                runner_hist_a + k * (runner_hist_b - runner_hist_a) / runner_hist_N ,
                runner_hist_a + (k + 1) * (runner_hist_b - runner_hist_a) / runner_hist_N ,
                (double)runner_hist_bins[k] );
    #endif
        
Pedro Gonnet's avatar
Pedro Gonnet committed
966
967
968
969
    /* Loop over the parts directly. */
    // for ( k = 0 ; k < N ; k++ )
    //     printf( " %i %e %e\n" , s.parts[k].id , s.parts[k].count , s.parts[k].count_dh );
    
970
    /* Dump the task data. */
971
972
973
974
975
976
977
    /* for ( k = 0 ; k < e.sched.nr_tasks ; k++ )
        if ( !e.sched.tasks[k].skip && !e.sched.tasks[k].implicit )
            printf( " %i %i %i %i %lli %lli %i %i\n" ,
                e.sched.tasks[k].rid , e.sched.tasks[k].type , e.sched.tasks[k].subtype , 
                (e.sched.tasks[k].cj == NULL) , e.sched.tasks[k].tic , e.sched.tasks[k].toc ,
                e.sched.tasks[k].ci->count , 
                (e.sched.tasks[k].cj==NULL)?0:e.sched.tasks[k].cj->count ); */
978
    
979
    /* Write final output. */
980
981
982
983
984
985
986
    if ( myrank == 0 )
        write_output( &e );
        
    #ifdef WITH_MPI
        if ( MPI_Finalize() != MPI_SUCCESS )
            error( "call to MPI_Finalize failed with error %i." , res );
    #endif
Pedro Gonnet's avatar
Pedro Gonnet committed
987
    
988
989
    /* Say goodbye. */
    message( "done." );
990
    
Pedro Gonnet's avatar
Pedro Gonnet committed
991
992
993
994
    /* All is calm, all is bright. */
    return 0;
    
    }