test.c 30.7 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

/* Ticks per second on this machine. */
#ifndef CPU_TPS
Pedro Gonnet's avatar
Pedro Gonnet committed
51
    #define CPU_TPS 2.40e9
Pedro Gonnet's avatar
Pedro Gonnet committed
52
53
#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
264
265

    }


/**
 * @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 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
266
267
268
269
270
    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
271
272
    
    /* Loop over all particle pairs. */
Pedro Gonnet's avatar
Pedro Gonnet committed
273
    #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
274
275
    for ( j = 0 ; j < N ; j++ ) {
        if ( j % 1000 == 0 ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
276
            printf( "pairs_n2: j=%i.\n" , j );
Pedro Gonnet's avatar
Pedro Gonnet committed
277
278
279
280
281
282
283
284
285
286
287
288
            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
289
            r2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2];
Pedro Gonnet's avatar
Pedro Gonnet committed
290
            if ( r2 < parts[j].h*parts[j].h || r2 < parts[k].h*parts[k].h ) {
291
                runner_iact_density( r2 , NULL , parts[j].h , parts[k].h , &parts[j] , &parts[k] );
Pedro Gonnet's avatar
Pedro Gonnet committed
292
                /* if ( parts[j].h / parts[k].h > maxratio )
Pedro Gonnet's avatar
Pedro Gonnet committed
293
294
                    #pragma omp critical
                    {
Pedro Gonnet's avatar
Pedro Gonnet committed
295
                    maxratio = parts[j].h / parts[k].h;
Pedro Gonnet's avatar
Pedro Gonnet committed
296
297
                    mj = j; mk = k;
                    }
Pedro Gonnet's avatar
Pedro Gonnet committed
298
                else if ( parts[k].h / parts[j].h > maxratio )
Pedro Gonnet's avatar
Pedro Gonnet committed
299
300
                    #pragma omp critical
                    {
Pedro Gonnet's avatar
Pedro Gonnet committed
301
                    maxratio = parts[k].h / parts[j].h;
Pedro Gonnet's avatar
Pedro Gonnet committed
302
                    mj = j; mk = k;
Pedro Gonnet's avatar
Pedro Gonnet committed
303
                    } */
Pedro Gonnet's avatar
Pedro Gonnet committed
304
305
306
                }
            }
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
307
308
309
        
    /* Aggregate the results. */
    for ( k = 0 ; k < N ; k++ ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
310
        // count += parts[k].icount;
311
312
313
        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
314
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
315
316
            
    /* Dump the result. */
Pedro Gonnet's avatar
Pedro Gonnet committed
317
    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
318
319
    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
320
321
        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
322
323
324
325
326
327
        parts[mj].h , parts[mk].h ); fflush(stdout); */
    fflush(stdout);
            
    }


Pedro Gonnet's avatar
Pedro Gonnet committed
328
void pairs_single_density ( double *dim , long long int pid , struct part *__restrict__ parts , int N , int periodic ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
329
330
331
332
333

    int i, k;
    // int mj, mk;
    // double maxratio = 1.0;
    double r2, dx[3];
334
335
336
    float fdx[3];
    struct part p;
    // double ih = 12.0/6.25;
337
338
339
340
341
    
    /* Find "our" part. */
    for ( k = 0 ; k < N && parts[k].id != pid ; k++ );
    if ( k == N )
        error( "Part not found." );
342
    p = parts[k];
343
    printf( "pairs_single: part[%i].id == %lli.\n" , k , pid );
Pedro Gonnet's avatar
Pedro Gonnet committed
344
    
345
346
347
348
349
    p.rho = 0.0;
    p.density.wcount = 0.0;
    // p.icount = 0;
    p.rho_dh = 0.0;
            
Pedro Gonnet's avatar
Pedro Gonnet committed
350
351
    /* Loop over all particle pairs. */
    for ( k = 0 ; k < N ; k++ ) {
352
        if ( parts[k].id == p.id )
Pedro Gonnet's avatar
Pedro Gonnet committed
353
354
            continue;
        for ( i = 0 ; i < 3 ; i++ ) {
355
            dx[i] = p.x[i] - parts[k].x[i];
Pedro Gonnet's avatar
Pedro Gonnet committed
356
357
358
359
360
361
            if ( periodic ) {
                if ( dx[i] < -dim[i]/2 )
                    dx[i] += dim[i];
                else if ( dx[i] > dim[i]/2 )
                    dx[i] -= dim[i];
                }
362
            fdx[i] = dx[i];
Pedro Gonnet's avatar
Pedro Gonnet committed
363
            }
364
365
366
367
368
        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) ,
369
                parts[k].id , (int)(parts[k].x[0]*ih) , (int)(parts[k].x[1]*ih) , (int)(parts[k].x[2]*ih) ,
370
                sqrtf(r2) ); */
Pedro Gonnet's avatar
Pedro Gonnet committed
371
372
373
374
            }
        }
        
    /* Dump the result. */
375
    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
376
    fflush(stdout);
Pedro Gonnet's avatar
Pedro Gonnet committed
377
    
Pedro Gonnet's avatar
Pedro Gonnet committed
378
379
380
    }


Pedro Gonnet's avatar
Pedro Gonnet committed
381
void pairs_single_grav ( double *dim , long long int pid , struct gpart *__restrict__ parts , int N , int periodic ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
382

Pedro Gonnet's avatar
Pedro Gonnet committed
383
384
385
386
387
388
389
    int i, k;
    // int mj, mk;
    // double maxratio = 1.0;
    double r2, dx[3];
    float fdx[3], a[3] = { 0.0 , 0.0 , 0.0 }, aabs[3] = { 0.0 , 0.0 , 0.0 };
    struct gpart pi, pj;
    // double ih = 12.0/6.25;
Pedro Gonnet's avatar
Pedro Gonnet committed
390
    
Pedro Gonnet's avatar
Pedro Gonnet committed
391
392
393
394
395
396
397
398
    /* Find "our" part. */
    for ( k = 0 ; k < N ; k++ )
        if ( ( parts[k].id > 0 && parts[k].part->id == pid ) || parts[k].id == -pid )
            break;
    if ( k == N )
        error( "Part not found." );
    pi = parts[k];
    pi.a[0] = 0.0f; pi.a[1] = 0.0f; pi.a[2] = 0.0f;
Pedro Gonnet's avatar
Pedro Gonnet committed
399
400
401
    
    /* Loop over all particle pairs. */
    for ( k = 0 ; k < N ; k++ ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
402
        if ( parts[k].id == pi.id )
Pedro Gonnet's avatar
Pedro Gonnet committed
403
            continue;
Pedro Gonnet's avatar
Pedro Gonnet committed
404
        pj = parts[k];
Pedro Gonnet's avatar
Pedro Gonnet committed
405
        for ( i = 0 ; i < 3 ; i++ ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
406
            dx[i] = pi.x[i] - pj.x[i];
Pedro Gonnet's avatar
Pedro Gonnet committed
407
408
409
410
411
412
            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
413
            fdx[i] = dx[i];
Pedro Gonnet's avatar
Pedro Gonnet committed
414
            }
Pedro Gonnet's avatar
Pedro Gonnet committed
415
416
417
418
419
        r2 = fdx[0]*fdx[0] + fdx[1]*fdx[1] + fdx[2]*fdx[2];
        runner_iact_grav( r2 , fdx , &pi , &pj );
        a[0] += pi.a[0]; a[1] += pi.a[1]; a[2] += pi.a[2];
        aabs[0] += fabsf( pi.a[0] ); aabs[1] += fabsf( pi.a[1] ); aabs[2] += fabsf( pi.a[2] );
        pi.a[0] = 0.0f; pi.a[1] = 0.0f; pi.a[2] = 0.0f;
Pedro Gonnet's avatar
Pedro Gonnet committed
420
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
421
422
423
        
    /* Dump the result. */
    message( "acceleration on gpart %lli is a=[ %e %e %e ], |a|=[ %.2e %.2e %.2e ].\n" , pi.part->id , a[0] , a[1] , a[2] , aabs[0] , aabs[1] , aabs[2] );
424
    
Pedro Gonnet's avatar
Pedro Gonnet committed
425
426
427
    }


428
429
430
431
432
433
434
435
436
437
/**
 * @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;
438
439
440
    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)));
441
442

    for ( k = 0 ; k <= N ; k++ ) {
443
        x = ((float)k) / N;
Pedro Gonnet's avatar
Pedro Gonnet committed
444
        x4[3] = x4[2]; x4[2] = x4[1]; x4[1] = x4[0]; x4[0] = x;
445
        kernel_deval( x , &w , &dw_dx );
446
        // kernel_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
Pedro Gonnet's avatar
Pedro Gonnet committed
447
448
449
450
451
452
        printf( " %e %e %e %e %e %e %e\n" , x , w , dw_dx , w4[0] , w4[1] , w4[2] , w4[3] );
        }

    }


Pedro Gonnet's avatar
Pedro Gonnet committed
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
void gravity_dump ( float r_max , int N ) {

    int k;
    float x, w;
    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)));
    
    float gadget ( float r ) {
        float fac, h_inv, u, r2 = r*r;
        if ( r >= const_epsilon )
            fac = 1.0f / (r2 * r);
        else {
            h_inv = 1. / const_epsilon;
            u = r * h_inv;
            if ( u < 0.5 )
                fac = const_iepsilon3 * (10.666666666667 + u * u * (32.0 * u - 38.4));
            else
                fac = const_iepsilon3 * (21.333333333333 - 48.0 * u +
                                     38.4 * u * u - 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
            }
        return const_G * fac;
        }

    for ( k = 1 ; k <= N ; k++ ) {
        x = (r_max * k) / N;
        x4[3] = x4[2]; x4[2] = x4[1]; x4[1] = x4[0]; x4[0] = x;
        kernel_grav_eval( x , &w );
        w *= const_G / ( x*x*x );
        // blender_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
        printf( " %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n" , x , w*x , w4[0] , w4[1] , w4[2] , w4[3] , gadget(x)*x );
        }

    }


Pedro Gonnet's avatar
Pedro Gonnet committed
489
490
491
492
493
494
495
496
497
498
499
500
501
502
/**
 * @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++ ) {
503
504
        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
505
506
507
508
509
510
511
512
513
        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;
514
        Pi[0].density.wcount = 0; Pj[0].density.wcount = 0;
Pedro Gonnet's avatar
Pedro Gonnet committed
515
        runner_iact_density( r2[0] , NULL , hi[0] , hj[0] , &Pi[0] , &Pj[0] );
516
517
518
519
520
        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
521
        runner_iact_vec_density( r2 , NULL , hi , hj , pi , pj );
522
        printf( " %e %e %e %e\n" , Pi[0].density.wcount , Pi[1].density.wcount , Pi[2].density.wcount , Pi[3].density.wcount );
523
524
525
        }

    }
Pedro Gonnet's avatar
Pedro Gonnet committed
526
527
528
529
530
531
532
533
534


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

Pedro Gonnet's avatar
Pedro Gonnet committed
535
    int c, icount, j, k, N = -1, periodic = 1;
Pedro Gonnet's avatar
Pedro Gonnet committed
536
    int nr_threads = 1, nr_queues = -1, runs = INT_MAX;
537
    int data[2];
Pedro Gonnet's avatar
Pedro Gonnet committed
538
    double dim[3] = { 1.0 , 1.0 , 1.0 }, shift[3] = { 0.0 , 0.0 , 0.0 };
539
    double h_max = -1.0 , scaling = 1.0;
Pedro Gonnet's avatar
Pedro Gonnet committed
540
    double clock = DBL_MAX;
541
    struct part *parts = NULL;
Pedro Gonnet's avatar
Pedro Gonnet committed
542
    struct space s;
Pedro Gonnet's avatar
Pedro Gonnet committed
543
    struct engine e;
544
    struct UnitSystem us;
545
    char ICfileName[200];
546
    float dt_max = 0.0f;
Pedro Gonnet's avatar
Pedro Gonnet committed
547
    ticks tic;
548
    int nr_nodes = 1, myrank = 0, grid[3] = { 1 , 1 , 1 };
Pedro Gonnet's avatar
Pedro Gonnet committed
549
    
550
    /* Choke on FP-exceptions. */
Pedro Gonnet's avatar
Pedro Gonnet committed
551
    // feenableexcept( FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW );
552
    
Pedro Gonnet's avatar
Pedro Gonnet committed
553
554
555
556
    /* Just dump the gravity potential and leave. */
    /* gravity_dump( 0.005 , 1000 );
    return 0; */
    
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
#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
574
575
576

    /* Greeting message */
    message( "This is SWIFT version %s\n", git_revision());
577
    
Pedro Gonnet's avatar
Pedro Gonnet committed
578
579
    /* Init the space. */
    bzero( &s , sizeof(struct space) );
580
581

    /* Parse the options */
582
    while ( ( c = getopt( argc , argv  , "a:c:d:f:g:m:q:r:s:t:w:z:" ) ) != -1 )
583
584
585
586
587
      switch( c )
	{
	case 'a':
	  if ( sscanf( optarg , "%lf" , &scaling ) != 1 )
	    error( "Error parsing cutoff scaling." );
588
589
	  if ( myrank == 0 )
        message( "scaling cutoff by %.3f." , scaling ); fflush(stdout);
590
	  break;
Pedro Gonnet's avatar
Pedro Gonnet committed
591
592
593
	case 'c':
	  if ( sscanf( optarg , "%lf" , &clock ) != 1 )
	    error( "Error parsing clock." );
594
595
	  if ( myrank == 0 )
        message( "clock set to %.3e." , clock ); fflush(stdout);
Pedro Gonnet's avatar
Pedro Gonnet committed
596
	  break;
597
	case 'd':
598
	  if ( sscanf( optarg , "%f" , &dt_max ) != 1 )
599
	    error( "Error parsing timestep." );
600
601
	  if ( myrank == 0 )
        message( "dt set to %e." , dt_max ); fflush(stdout);
602
603
604
605
	  break;
	case 'f':
	  if( !strcpy(ICfileName, optarg))
	    error("Error parsing IC file name.");
Pedro Gonnet's avatar
Pedro Gonnet committed
606
607
	  // if ( myrank == 0 )
        // message("IC to be read from file '%s'.", ICfileName);
608
609
610
611
612
613
	  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);
614
	  break;
615
616
617
	case 'm':
	  if ( sscanf( optarg , "%lf" , &h_max ) != 1 )
	    error( "Error parsing h_max." );
618
619
	  if ( myrank == 0 )
        message( "maximum h set to %e." , h_max ); fflush(stdout);
620
621
622
623
624
625
626
627
628
629
630
631
	  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." );
632
633
	  if ( myrank == 0 )
        message( "will shift parts by [ %.3f %.3f %.3f ]." , shift[0] , shift[1] , shift[2] );
634
635
636
637
638
639
	  break;
	case 't':
	  if ( sscanf( optarg , "%d" , &nr_threads ) != 1 )
	    error( "Error parsing number of threads." );
	  omp_set_num_threads( nr_threads );
	  break;
640
641
642
	case 'w':
	  if ( sscanf( optarg , "%d" , &space_subsize ) != 1 )
	    error( "Error parsing sub size." );
643
644
	  if ( myrank == 0 )
        message( "sub size set to %i." , space_subsize );
645
	  break;
646
647
648
	case 'z':
	  if ( sscanf( optarg , "%d" , &space_splitsize ) != 1 )
	    error( "Error parsing split size." );
649
650
	  if ( myrank == 0 )
        message( "split size set to %i." , space_splitsize );
651
652
653
654
655
656
	  break;
	case '?':
	  error( "Unknown option." );
	  break;

	}
Pedro Gonnet's avatar
Pedro Gonnet committed
657
    
658

659
    /* How large are the parts? */
Pedro Gonnet's avatar
Pedro Gonnet committed
660
    if ( myrank == 0 ) {
661
        message( "sizeof(struct part) is %li bytes." , (long int)sizeof( struct part ));
Pedro Gonnet's avatar
Pedro Gonnet committed
662
663
        message( "sizeof(struct gpart) is %li bytes." , (long int)sizeof( struct gpart ));
        }
664

665
666
667
668
669
    /* Initilaize unit system */
    initUnitSystem(&us);
    if ( myrank == 0 )
      {
	message( "Unit system: U_M = %e g.", us.UnitMass_in_cgs );
670
	message( "Unit system: U_L = %e cm.", us.UnitLength_in_cgs );
671
672
673
	message( "Unit system: U_t = %e s.", us.UnitTime_in_cgs );
	message( "Unit system: U_I = %e A.", us.UnitCurrent_in_cgs );
	message( "Unit system: U_T = %e K.", us.UnitTemperature_in_cgs );
674
675
	message( "Density units: %e a^%f h^%f.", conversionFactor(&us, UNIT_CONV_DENSITY), aFactor(&us, UNIT_CONV_DENSITY), hFactor(&us, UNIT_CONV_DENSITY) );
	message( "Entropy units: %e a^%f h^%f.", conversionFactor(&us, UNIT_CONV_ENTROPY), aFactor(&us, UNIT_CONV_ENTROPY), hFactor(&us, UNIT_CONV_ENTROPY) );
676
677
      }

678
679
    /* Read particles and space information from (GADGET) IC */
    tic = getticks();
680
#ifdef WITH_MPI
681
    read_ic_parallel( ICfileName , dim , &parts , &N , &periodic, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL );
682
#else
683
    read_ic( ICfileName , dim , &parts , &N , &periodic );
684
#endif
685
686
    if ( myrank == 0 )
        message( "reading particle properties took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
687
688
    
    /* Apply h scaling */
Pedro Gonnet's avatar
Pedro Gonnet committed
689
    if( scaling != 1.0 )
690
      for ( k = 0 ; k < N ; k++ )
691
	    parts[k].h *= scaling;
692

693
694
695
    /* 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
696
697
698
	    parts[k].x[0] += shift[0];
	    parts[k].x[1] += shift[1];
	    parts[k].x[2] += shift[2];
699
700
      }

Pedro Gonnet's avatar
Pedro Gonnet committed
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
    /* printParticle( parts , 10312237508790 , N );
    printParticle( parts , 10312286091950 , N ); */
    /* for ( k = 0 ; k < N ; k++ )
	if ( parts[k].id == 10312286091950 ||
	     parts[k].id == 10286889371446  ||
	     parts[k].id == 9536045071298 ||
	     parts[k].id == 12726778692106  ||
	     parts[k].id == 9479892852626  ||
	     parts[k].id == 9535843125514  ||
	     parts[k].id == 14151507889834  ||
	     parts[k].id == 14144038209438  ||
	     parts[k].id == 14121890205050  ||
	     parts[k].id == 5868762382714  ||
	     parts[k].id == 12840527117206 ||
	     parts[k].id == 10292087642778  ||
	     parts[k].id == 9465178320650  ||
	     parts[k].id == 2834846537770  ||
	     parts[k].id == 9483000048314  ||
	     parts[k].id == 10247332828902  ||
	     parts[k].id == 10223834653674  ||
	     parts[k].id == 16719632108962  ||
	     parts[k].id == 16759192850622  ||
	     parts[k].id == 9483599082554  ||
	     parts[k].id == 10247340329226 )
            parts[k] = parts[--N]; */

Pedro Gonnet's avatar
Pedro Gonnet committed
727
728
729
730
731
732
733
    /* 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 );
734
735
    if ( myrank == 0 )
        message( "space_init took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
736

737
    /* Set the default time step to 1.0f. */
738
739
    if ( myrank == 0 )
        message( "dt_max is %e." , dt_max );
740
    
Pedro Gonnet's avatar
Pedro Gonnet committed
741
    /* Say a few nice things about the space we just created. */
742
743
744
745
746
747
748
749
    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
750
    
Pedro Gonnet's avatar
Pedro Gonnet committed
751
    /* Verify that each particle is in it's propper cell. */
752
753
754
755
756
    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
757
    
758
759
760
761
762
    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] );
        }
763
    
Pedro Gonnet's avatar
Pedro Gonnet committed
764
765
766
    /* Dump the particle positions. */
    // space_map_parts( &s , &map_dump , shift );
    
767
768
    
    /* Initialize the engine with this space. */
Pedro Gonnet's avatar
Pedro Gonnet committed
769
    tic = getticks();
Pedro Gonnet's avatar
Pedro Gonnet committed
770
    message( "nr_nodes is %i." , nr_nodes );
771
772
773
774
775
776
777
    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. */
    engine_split( &e , grid );
778
    engine_redistribute ( &e );
779
#endif
780
781

    /* Write the state of the system as it is before starting time integration. */
782
783
    tic = getticks();
#ifdef WITH_MPI
784
    write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
785
#else
786
    write_output(&e, &us);
787
788
#endif
    message( "writing particle properties took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
Pedro Gonnet's avatar
Pedro Gonnet committed
789
    
790
791
792
793
794
795
    /* 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
796
797
    /* Inauguration speech. */
    if ( runs < INT_MAX )
798
        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
799
    else
800
        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
801
802
    fflush(stdout);
    
Pedro Gonnet's avatar
Pedro Gonnet committed
803
804
805
806
807
808
809
810
    /* Set a target particle. */
    /* long long int pid[5];
    unsigned int seed = 6178;
    for ( k = 0 ; k < 5 ; k++ )
        pid[k] = s.gparts[ rand_r( &seed ) % N ].part->id;
    for ( k = 0 ; k < 5 ; k++ )
        pairs_single_grav( dim , pid[k] , s.gparts , N , 0 ); */
    
Pedro Gonnet's avatar
Pedro Gonnet committed
811
    /* Legend. */
812
813
    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
814
    
Pedro Gonnet's avatar
Pedro Gonnet committed
815
    /* Let loose a runner on the space. */
Pedro Gonnet's avatar
Pedro Gonnet committed
816
    for ( j = 0 ; j < runs && e.time < clock ; j++ ) {
817
    
818
819
        /* Repartition the space amongst the nodes? */
        #if defined(WITH_MPI) && defined(HAVE_METIS)
820
            if ( j == 2 )
821
822
823
                e.forcerepart = 1;
        #endif
        
824
825
826
827
        /* Force a rebuild for testing. */
        /* if ( j % 4 == 1 )
            e.forcerebuild = 1; */
        
828
        // 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);
829
        timers_reset( timers_mask_all );
Pedro Gonnet's avatar
Pedro Gonnet committed
830
831
832
833
        #ifdef COUNTER
            for ( k = 0 ; k < runner_counter_count ; k++ )
                runner_counter[k] = 0;
        #endif
834
835
        
        /* Take a step. */
Pedro Gonnet's avatar
Pedro Gonnet committed
836
        engine_step( &e );
Pedro Gonnet's avatar
Pedro Gonnet committed
837
        
838
839
840
        if ( j % 100 == 0 )
	  {
#ifdef WITH_MPI
841
             write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
842
#else
843
             write_output(&e, &us);
844
845
#endif
          }
846
                
Pedro Gonnet's avatar
Pedro Gonnet committed
847
        /* Dump a line of agregate output. */
848
849
850
851
852
853
854
855
856
857
        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);
            }
Pedro Gonnet's avatar
Pedro Gonnet committed
858
859
860
        /* for ( k = 0 ; k < 5 ; k++ )
            printgParticle( s.gparts , pid[k] , N ); */
        
Pedro Gonnet's avatar
Pedro Gonnet committed
861
862
        }
        
863
864
865
866
867
868
869
870
871
    /* 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
872
873

    // write_output( &e );
874
        
Pedro Gonnet's avatar
Pedro Gonnet committed
875
876
877
878
    /* 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 );
    
879
    /* Dump the task data. */
Pedro Gonnet's avatar
Pedro Gonnet committed
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
    /* #ifdef WITH_MPI
    for ( j = 0 ; j < nr_nodes ; j++ ) {
    MPI_Barrier( MPI_COMM_WORLD );
    if ( j == myrank ) {
    printf( " %03i 0 0 0 0 %lli 0 0 0 0\n" , myrank , e.tic_step );
    for ( k = 0 ; k < e.sched.nr_tasks ; k++ )
        if ( !e.sched.tasks[k].skip && !e.sched.tasks[k].implicit )
            printf( " %03i %i %i %i %i %lli %lli %i %i %i\n" ,
                myrank ,
                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)?e.sched.tasks[k].cj->count:0 , e.sched.tasks[k].flags); 
    fflush(stdout);
    sleep(1);
    }
    }
    #else
    for ( k = 0 ; k < e.sched.nr_tasks ; k++ )
898
        if ( !e.sched.tasks[k].skip && !e.sched.tasks[k].implicit )
Pedro Gonnet's avatar
Pedro Gonnet committed
899
                printf( " %i %i %i %i %lli %lli %i %i\n" ,
900
901
902
                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 , 
Pedro Gonnet's avatar
Pedro Gonnet committed
903
904
                (e.sched.tasks[k].cj==NULL)?0:e.sched.tasks[k].cj->count ); 
    #endif */
905
    
906
    /* Write final output. */
907
#ifdef WITH_MPI
908
	write_output_parallel( &e, &us, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL );
909
#else
910
	write_output( &e , &us );
911
#endif
912
        
913
#ifdef WITH_MPI
914
915
916
        if ( MPI_Finalize() != MPI_SUCCESS )
            error( "call to MPI_Finalize failed with error %i." , res );
    #endif
Pedro Gonnet's avatar
Pedro Gonnet committed
917
    
918
919
    /* Say goodbye. */
    message( "done." );
920
    
Pedro Gonnet's avatar
Pedro Gonnet committed
921
922
923
924
    /* All is calm, all is bright. */
    return 0;
    
    }