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

    /* Greeting message */
576
    message( "This is %s\n", package_description() );
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
681
#if defined( WITH_MPI ) 
#if defined( HAVE_PARALLEL_HDF5 )
682
    read_ic_parallel( ICfileName , dim , &parts , &N , &periodic, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL );
683
#else
684
685
686
687
    read_ic_serial( ICfileName , dim , &parts , &N , &periodic, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL );
#endif
#else
    read_ic_single( ICfileName , dim , &parts , &N , &periodic );
688
#endif
689

690
691
    if ( myrank == 0 )
        message( "reading particle properties took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
692
693
    
    /* Apply h scaling */
Pedro Gonnet's avatar
Pedro Gonnet committed
694
    if( scaling != 1.0 )
695
      for ( k = 0 ; k < N ; k++ )
696
	    parts[k].h *= scaling;
697

698
699
700
    /* 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
701
702
703
	    parts[k].x[0] += shift[0];
	    parts[k].x[1] += shift[1];
	    parts[k].x[2] += shift[2];
704
705
      }

Pedro Gonnet's avatar
Pedro Gonnet committed
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
    /* 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
732
733
734
735
736
737
738
    /* 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 );
739
740
    if ( myrank == 0 )
        message( "space_init took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
741

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

786
    message("Before write !");
787
    /* Write the state of the system as it is before starting time integration. */
788
    tic = getticks();
789
790
#if defined( WITH_MPI ) 
#if defined( HAVE_PARALLEL_HDF5 )
791
    write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
792
#else
793
794
795
796
    write_output_serial(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
#endif
#else
    write_output_single(&e, &us);
797
798
#endif
    message( "writing particle properties took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
Pedro Gonnet's avatar
Pedro Gonnet committed
799
    
800
801
802
803
804
805
    /* 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
806
807
    /* Inauguration speech. */
    if ( runs < INT_MAX )
808
        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
809
    else
810
        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
811
812
    fflush(stdout);
    
Pedro Gonnet's avatar
Pedro Gonnet committed
813
814
815
816
817
818
819
820
    /* 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
821
    /* Legend. */
822
823
    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
824
    
Pedro Gonnet's avatar
Pedro Gonnet committed
825
    /* Let loose a runner on the space. */
Pedro Gonnet's avatar
Pedro Gonnet committed
826
    for ( j = 0 ; j < runs && e.time < clock ; j++ ) {
827
    
828
829
        /* Repartition the space amongst the nodes? */
        #if defined(WITH_MPI) && defined(HAVE_METIS)
830
            if ( j % 100 == 2 )
831
832
833
                e.forcerepart = 1;
        #endif
        
834
        /* Force a rebuild for testing. */
835
        /* if ( j % 4 == 3 )
836
837
            e.forcerebuild = 1; */
        
838
        // 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);
839
        timers_reset( timers_mask_all );
Pedro Gonnet's avatar
Pedro Gonnet committed
840
841
842
843
        #ifdef COUNTER
            for ( k = 0 ; k < runner_counter_count ; k++ )
                runner_counter[k] = 0;
        #endif
844
845
        
        /* Take a step. */
Pedro Gonnet's avatar
Pedro Gonnet committed
846
        engine_step( &e );
Pedro Gonnet's avatar
Pedro Gonnet committed
847
        
848
849
        if ( j % 100 == 0 )
	  {
850
851
852
853
854
855
856

#if defined( WITH_MPI ) 
#if defined( HAVE_PARALLEL_HDF5 )
	    write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
#else
	    write_output_serial(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
#endif
857
#else
858
	    write_output_single(&e, &us);
859
#endif
860

861
          }
862
                
Pedro Gonnet's avatar
Pedro Gonnet committed
863
        /* Dump a line of agregate output. */
864
865
866
867
868
869
870
871
872
873
        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
874
875
876
        /* for ( k = 0 ; k < 5 ; k++ )
            printgParticle( s.gparts , pid[k] , N ); */
        
877
    }
Pedro Gonnet's avatar
Pedro Gonnet committed
878
        
879
    /* Print the values of the runner histogram. */
880
881
882
883
884
885
886
887
#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
888

Pedro Gonnet's avatar
Pedro Gonnet committed
889
890
891
892
    /* 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 );
    
893
    /* Dump the task data. */
Pedro Gonnet's avatar
Pedro Gonnet committed
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
    /* #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++ )
912
        if ( !e.sched.tasks[k].skip && !e.sched.tasks[k].implicit )
Pedro Gonnet's avatar
Pedro Gonnet committed
913
                printf( " %i %i %i %i %lli %lli %i %i\n" ,
914
915
916
                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
917
918
                (e.sched.tasks[k].cj==NULL)?0:e.sched.tasks[k].cj->count ); 
    #endif */
919
    
920
    /* Write final output. */
921
922
923
924
925
926
#if defined( WITH_MPI ) 
#if defined( HAVE_PARALLEL_HDF5 )
    write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
#else
    write_output_serial(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
#endif
927
#else
928
    write_output_single(&e, &us);
929
#endif
930
        
931
#ifdef WITH_MPI
932
933
934
    if ( MPI_Finalize() != MPI_SUCCESS )
      error( "call to MPI_Finalize failed with error %i." , res );
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
935
    
936
937
    /* Say goodbye. */
    message( "done." );
938
    
Pedro Gonnet's avatar
Pedro Gonnet committed
939
940
941
    /* All is calm, all is bright. */
    return 0;
    
942
}