engine.c 77.5 KB
Newer Older
Pedro Gonnet's avatar
Pedro Gonnet committed
1
/*******************************************************************************
2
 * This file is part of SWIFT.
Pedro Gonnet's avatar
Pedro Gonnet committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
 * 
 * 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/>.
 * 
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
#include <stdio.h>
#include <stdlib.h>
26
#include <unistd.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
27
28
29
30
31
32
33
34
#include <string.h>
#include <pthread.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include <omp.h>
#include <sched.h>

35
36
37
38
39
/* MPI headers. */
#ifdef WITH_MPI
    #include <mpi.h>
#endif

40
41
42
43
44
/* METIS headers. */
#ifdef HAVE_METIS
    #include <metis.h>
#endif

Pedro Gonnet's avatar
Pedro Gonnet committed
45
/* Local headers. */
46
#include "const.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
47
#include "cycle.h"
48
#include "atomic.h"
49
#include "timers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
50
#include "const.h"
51
#include "vector.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
52
53
54
#include "lock.h"
#include "task.h"
#include "part.h"
55
#include "debug.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
56
#include "space.h"
57
#include "multipole.h"
58
#include "cell.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
59
#include "queue.h"
60
#include "scheduler.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
61
62
#include "engine.h"
#include "runner.h"
63
#include "proxy.h"
64
#include "error.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
65

66
67
68
69
70
71
72
#ifdef LEGACY_GADGET2_SPH
#include "runner_iact_legacy.h"
#else
#include "runner_iact.h"
#endif


Pedro Gonnet's avatar
Pedro Gonnet committed
73
74
75
76
/* Convert cell location to ID. */
#define cell_getid( cdim , i , j , k ) ( (int)(k) + (cdim)[2]*( (int)(j) + (cdim)[1]*(int)(i) ) )


77
78
79
80
/** The rank of the engine as a global variable (for messages). */
int engine_rank;


81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
/**
 * @brief Link a density/force task to a cell.
 *
 * @param e The #engine.
 * @param l The #link.
 * @param t The #task.
 *
 * @return The new #link pointer.
 */
 
struct link *engine_addlink( struct engine *e , struct link *l , struct task *t ) {

    struct link *res = &e->links[ atomic_inc( &e->nr_links ) ];
    res->next = l;
    res->t = t;
    return res;

    }


/**
 * @brief Generate the ghost and kick tasks for a hierarchy of cells.
 *
 * @param e The #engine.
 * @param c The #cell.
 * @param super The super #cell.
 */
 
void engine_mkghosts ( struct engine *e , struct cell *c , struct cell *super ) {

    int k;
    struct scheduler *s = &e->sched;

    /* Am I the super-cell? */
    if ( super == NULL && c->nr_tasks > 0 ) {
    
        /* Remember me. */
        super = c;
        
        /* Local tasks only... */
        if ( c->nodeID == e->nodeID ) {
        
            /* Generate the ghost task. */
            c->ghost = scheduler_addtask( s , task_type_ghost , task_subtype_none , 0 , 0 , c , NULL , 0 );

            /* Add the kick2 task. */
            c->kick2 = scheduler_addtask( s , task_type_kick2 , task_subtype_none , 0 , 0 , c , NULL , 0 );

            /* Add the kick1 task if needed. */
            if ( !(e->policy & engine_policy_fixdt) )
                c->kick1 = scheduler_addtask( s , task_type_kick1 , task_subtype_none , 0 , 0 , c , NULL , 0 );
                
            }
            
        }
        
    /* Set the super-cell. */
    c->super = super;
        
    /* Recurse. */
    if ( c->split )
        for ( k = 0 ; k < 8 ; k++ )
            if ( c->progeny[k] != NULL )
                engine_mkghosts( e , c->progeny[k] , super );
    
    }


149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
/**
 * @brief Redistribute the particles amongst the nodes accorind
 *      to their cell's node IDs.
 *
 * @param e The #engine.
 */
 
void engine_redistribute ( struct engine *e ) {

#ifdef WITH_MPI

    int i, j, k, cid;
    int nr_nodes = e->nr_nodes, nodeID = e->nodeID;
    struct space *s = e->s;
    int my_cells = 0;
    int *cdim = s->cdim;
    struct cell *cells = s->cells;
166
    int nr_cells = s->nr_cells;
167
168
169
170
171

    /* Start by sorting the particles according to their nodes and
       getting the counts. */
    int *counts, *dest;
    struct part *parts = s->parts;
Pedro Gonnet's avatar
Pedro Gonnet committed
172
    double ih[3];
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
    ih[0] = s->ih[0]; ih[1] = s->ih[1]; ih[2] = s->ih[2];
    if ( ( counts = (int *)malloc( sizeof(int) * nr_nodes * nr_nodes ) ) == NULL ||
         ( dest = (int *)malloc( sizeof(int) * s->nr_parts ) ) == NULL )
        error( "Failed to allocate count and dest buffers." );
    bzero( counts , sizeof(int) * nr_nodes * nr_nodes );
    for ( k = 0 ; k < s->nr_parts ; k++ ) {
        cid = cell_getid( cdim , parts[k].x[0]*ih[0] , parts[k].x[1]*ih[1] , parts[k].x[2]*ih[2] );
        dest[k] = cells[ cid ].nodeID;
        counts[ nodeID*nr_nodes + dest[k] ] += 1;
        }
    parts_sort( s->parts , s->xparts , dest , s->nr_parts , 0 , nr_nodes-1 );
    
    /* Get all the counts from all the nodes. */
    if ( MPI_Allreduce( MPI_IN_PLACE , counts , nr_nodes * nr_nodes , MPI_INT , MPI_SUM , MPI_COMM_WORLD ) != MPI_SUCCESS )
        error( "Failed to allreduce particle transfer counts." );
188

189
190
191
192
193
194
    /* Get the new number of parts for this node, be generous in allocating. */
    int nr_parts = 0;
    for ( k = 0 ; k < nr_nodes ; k++ )
        nr_parts += counts[ k*nr_nodes + nodeID ];
    struct part *parts_new;
    struct xpart *xparts_new, *xparts = s->xparts;
195
196
    if ( posix_memalign( (void **)&parts_new , part_align , sizeof(struct part) * nr_parts * 1.2 ) != 0 ||
         posix_memalign( (void **)&xparts_new , part_align , sizeof(struct xpart) * nr_parts * 1.2 ) != 0 )
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
        error( "Failed to allocate new part data." );
        
    /* Emit the sends and recvs for the particle data. */
    MPI_Request *reqs;
    if ( ( reqs = (MPI_Request *)malloc( sizeof(MPI_Request) * 4 * nr_nodes ) ) == NULL )
        error( "Failed to allocate MPI request list." );
    for ( k = 0 ; k < 4*nr_nodes ; k++ )
        reqs[k] = MPI_REQUEST_NULL;
    for ( i = 0 , j = 0 , k = 0 ; k < nr_nodes ; k++ ) {
        if ( k == nodeID && counts[ nodeID*nr_nodes + k ] > 0 ) {
            memcpy( &parts_new[j] , &parts[i] , sizeof(struct part) * counts[ k*nr_nodes + nodeID ] );
            memcpy( &xparts_new[j] , &xparts[i] , sizeof(struct xpart) * counts[ k*nr_nodes + nodeID ] );
            i += counts[ nodeID*nr_nodes + k ];
            j += counts[ k*nr_nodes + nodeID ];
            }
        if ( k != nodeID && counts[ nodeID*nr_nodes + k ] > 0 ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
213
            if ( MPI_Isend( &parts[i] , sizeof(struct part) * counts[ nodeID*nr_nodes + k ] , MPI_BYTE , k , 2*(nodeID*nr_nodes + k) + 0 , MPI_COMM_WORLD , &reqs[4*k] ) != MPI_SUCCESS )
214
                error( "Failed to isend parts to node %i." , k );
Pedro Gonnet's avatar
Pedro Gonnet committed
215
            if ( MPI_Isend( &xparts[i] , sizeof(struct xpart) * counts[ nodeID*nr_nodes + k ] , MPI_BYTE , k , 2*(nodeID*nr_nodes + k) + 1 , MPI_COMM_WORLD , &reqs[4*k+1] ) != MPI_SUCCESS )
216
217
218
219
                error( "Failed to isend xparts to node %i." , k );
            i += counts[ nodeID*nr_nodes + k ];
            }
        if ( k != nodeID && counts[ k*nr_nodes + nodeID ] > 0 ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
220
            if ( MPI_Irecv( &parts_new[j] , sizeof(struct part) * counts[ k*nr_nodes + nodeID ] , MPI_BYTE , k , 2*(k*nr_nodes + nodeID) + 0 , MPI_COMM_WORLD , &reqs[4*k+2] ) != MPI_SUCCESS )
221
                error( "Failed to emit irecv of parts from node %i." , k );
Pedro Gonnet's avatar
Pedro Gonnet committed
222
            if ( MPI_Irecv( &xparts_new[j] , sizeof(struct xpart) * counts[ k*nr_nodes + nodeID ] , MPI_BYTE , k , 2*(k*nr_nodes + nodeID) + 1 , MPI_COMM_WORLD , &reqs[4*k+3] ) != MPI_SUCCESS )
223
224
225
226
227
                error( "Failed to emit irecv of parts from node %i." , k );
            j += counts[ k*nr_nodes + nodeID ];
            }
        }
        
Pedro Gonnet's avatar
Pedro Gonnet committed
228
    /* Wait for all the sends and recvs to tumble in. */
229
230
231
232
    MPI_Status stats[4*nr_nodes];
    int res;
    if ( ( res = MPI_Waitall( 4*nr_nodes , reqs , stats ) ) != MPI_SUCCESS ) {
        for ( k = 0 ; k < 4*nr_nodes ; k++ ) {
233
234
235
236
237
238
        char buff[ MPI_MAX_ERROR_STRING ];
        int res;
        MPI_Error_string( stats[k].MPI_ERROR , buff , &res );
        message( "request %i has error '%s'." , k , buff );
        }
    message( "counts is [ %i %i %i %i ]." , counts[0] , counts[1] , counts[2] , counts[3] );
239
        error( "Failed during waitall for part data." );
240
241
        }

242
243
244
245
246
247
248
249
250
251
252
253
254
    /* Verify that all parts are in the right place. */
    /* for ( k = 0 ; k < nr_parts ; k++ ) {
        cid = cell_getid( cdim , parts_new[k].x[0]*ih[0] , parts_new[k].x[1]*ih[1] , parts_new[k].x[2]*ih[2] );
        if ( cells[ cid ].nodeID != nodeID )
            error( "Received particle (%i) that does not belong here (nodeID=%i)." , k , cells[ cid ].nodeID );
        } */
        
    /* Set the new part data, free the old. */
    free( parts );
    free( xparts );
    s->parts = parts_new;
    s->xparts = xparts_new;
    s->nr_parts = nr_parts;
255
    s->size_parts = 1.2*nr_parts;
256
257
    
    /* Be verbose about what just happened. */
258
259
260
    for ( k = 0 ; k < nr_cells ; k++ )
        if ( cells[k].nodeID == nodeID )
            my_cells += 1;
261
262
263
264
265
266
267
268
269
270
271
272
273
274
    message( "node %i now has %i parts in %i cells." , nodeID , nr_parts , my_cells );
    
    /* Clean up other stuff. */
    free( reqs );
    free( counts );
    free( dest );
        
#else
    error( "SWIFT was not compiled with MPI and METIS support." );
#endif

    }


275
/**
276
 * @brief Repartition the cells amongst the nodes.
277
278
279
280
281
282
283
284
 *
 * @param e The #engine.
 */
 
void engine_repartition ( struct engine *e ) {

#if defined(WITH_MPI) && defined(HAVE_METIS)

285
    int i, j, k, l, cid, cjd, ii, jj, kk, res, w;
286
    idx_t *inds, *nodeIDs;
287
288
    idx_t *weights_v, *weights_e;
    struct space *s = e->s;
289
    int nr_cells = s->nr_cells, my_cells = 0;
290
291
292
293
    struct cell *cells = s->cells;
    int ind[3], *cdim = s->cdim;
    struct task *t, *tasks = e->sched.tasks;
    struct cell *ci, *cj;
294
    int nr_nodes = e->nr_nodes, nodeID = e->nodeID;
295
    float wscale = 0.0001, vscale = 0.001;
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
    
    /* Clear the repartition flag. */
    e->forcerepart = 0;
    
    /* Allocate the inds and weights. */
    if ( ( inds = (idx_t *)malloc( sizeof(idx_t) * 26*nr_cells ) ) == NULL ||
         ( weights_v = (idx_t *)malloc( sizeof(idx_t) * nr_cells ) ) == NULL ||
         ( weights_e = (idx_t *)malloc( sizeof(idx_t) * 26*nr_cells ) ) == NULL ||
         ( nodeIDs = (idx_t *)malloc( sizeof(idx_t) * nr_cells ) ) == NULL )
        error( "Failed to allocate inds and weights arrays." );
        
    /* Fill the inds array. */
    for ( cid = 0 ; cid < nr_cells ; cid++ ) {
        ind[0] = cells[cid].loc[0] / s->cells[cid].h[0] + 0.5;
        ind[1] = cells[cid].loc[1] / s->cells[cid].h[1] + 0.5;
        ind[2] = cells[cid].loc[2] / s->cells[cid].h[2] + 0.5;
        l = 0;
        for ( i = -1 ; i <= 1 ; i++ ) {
            ii = ind[0] + i;
            if ( ii < 0 ) ii += cdim[0];
            else if ( ii >= cdim[0] ) ii -= cdim[0];
            for ( j = -1 ; j <= 1 ; j++ ) {
                jj = ind[1] + j;
                if ( jj < 0 ) jj += cdim[1];
                else if ( jj >= cdim[1] ) jj -= cdim[1];
                for ( k = -1 ; k <= 1 ; k++ ) {
                    kk = ind[2] + k;
                    if ( kk < 0 ) kk += cdim[2];
                    else if ( kk >= cdim[2] ) kk -= cdim[2];
                    if ( i || j || k ) {
                        inds[ cid*26 + l ] = cell_getid( cdim , ii , jj , kk );
                        l += 1;
                        }
                    }
                }
            }
        }
        
    /* Init the weights arrays. */
335
336
    /* bzero( weights_e , sizeof(idx_t) * 26*nr_cells );
    bzero( weights_v , sizeof(idx_t) * nr_cells ); */
Pedro Gonnet's avatar
typo.    
Pedro Gonnet committed
337
    for ( k = 0 ; k < 26*nr_cells ; k++ )
338
        weights_e[k] = 1;
Pedro Gonnet's avatar
typo.    
Pedro Gonnet committed
339
    for ( k = 0 ; k < nr_cells ; k++ )
340
        weights_v[k] = 1;
341
342
343
344
345
346
347
348
349
350
    
    /* Loop over the tasks... */
    for ( j = 0 ; j < e->sched.nr_tasks ; j++ ) {
    
        /* Get a pointer to the kth task. */
        t = &tasks[j];
        
        /* Skip un-interesting tasks. */
        if ( t->type != task_type_self &&
             t->type != task_type_pair &&
351
352
353
354
             t->type != task_type_sub &&
             t->type != task_type_ghost &&
             t->type != task_type_kick1 &&
             t->type != task_type_kick2 )
355
            continue;
356
357
358
359
360
            
        /* Get the task weight. */
        w = ( t->toc - t->tic ) * wscale;
        if ( w < 0 )
            error( "Bad task weight (%i)." , w );
361
362
363
364
365
366
367
368
369
370
371
372
        
        /* Get the top-level cells involved. */
        for ( ci = t->ci ; ci->parent != NULL ; ci = ci->parent );
        if ( t->cj != NULL )
            for ( cj = t->cj ; cj->parent != NULL ; cj = cj->parent );
        else
            cj = NULL;
            
        /* Get the cell IDs. */
        cid = ci - cells;
            
        /* Different weights for different tasks. */
373
374
375
376
377
378
379
380
381
382
383
        if ( t->type == task_type_ghost ||
             t->type == task_type_kick1 || 
             t->type == task_type_kick2 ) {
             
            /* Particle updates add only to vertex weight. */
            weights_v[cid] += w;
            
            }
        
        /* Self interaction? */     
        else if ( ( t->type == task_type_self && ci->nodeID == nodeID ) ||
384
385
386
             ( t->type == task_type_sub && cj == NULL && ci->nodeID == nodeID ) ) {
        
            /* Self interactions add only to vertex weight. */
387
            weights_v[cid] += w;
388
389
390
391
392
393
394
395
396
397
398
            
            }
            
        /* Pair? */
        else if ( t->type == task_type_pair ||
                  ( t->type == task_type_sub && cj != NULL ) ) {
                  
            /* In-cell pair? */
            if ( ci == cj ) {
            
                /* Add weight to vertex for ci. */
399
                weights_v[cid] += w;
400
401
402
403
404
405
406
407
408
409
410
            
                }
                
            /* Distinct cells with local ci? */
            else if ( ci->nodeID == nodeID ) {
            
                /* Index of the jth cell. */
                cjd = cj - cells;
                
                /* Add half of weight to each cell. */
                if ( ci->nodeID == nodeID )
411
                    weights_v[cid] += 0.5 * w;
412
                if ( cj->nodeID == nodeID )
413
                    weights_v[cjd] += 0.5 * w;
414
415
416
                    
                /* Add Weight to edge. */
                for ( k = 26*cid ; inds[k] != cjd ; k++ );
417
                weights_e[ k ] += w;
418
                for ( k = 26*cjd ; inds[k] != cid ; k++ );
419
                weights_e[ k ] += w;
420
421
422
423
424
425
426
427
            
                }
                  
            }
    
        }
        
    /* Merge the weights arrays accross all nodes. */
428
#if IDXTYPEWIDTH==32
429
    if ( ( res = MPI_Reduce( ( nodeID == 0 ) ? MPI_IN_PLACE : weights_v , weights_v , nr_cells , MPI_INT , MPI_SUM , 0 , MPI_COMM_WORLD ) ) != MPI_SUCCESS ) {
430
431
432
#else
    if ( ( res = MPI_Reduce( ( nodeID == 0 ) ? MPI_IN_PLACE : weights_v , weights_v , nr_cells , MPI_LONG_LONG_INT , MPI_SUM , 0 , MPI_COMM_WORLD ) ) != MPI_SUCCESS ) {
#endif
433
434
435
436
        char buff[ MPI_MAX_ERROR_STRING ];
        MPI_Error_string( res , buff , &i );
        error( "Failed to allreduce vertex weights (%s)." , buff );
        }
437
#if IDXTYPEWIDTH==32
438
    if ( MPI_Reduce( ( nodeID == 0 ) ? MPI_IN_PLACE : weights_e , weights_e , 26*nr_cells , MPI_INT , MPI_SUM , 0 , MPI_COMM_WORLD ) != MPI_SUCCESS )
439
440
441
442
#else
    if ( MPI_Reduce( ( nodeID == 0 ) ? MPI_IN_PLACE : weights_e , weights_e , 26*nr_cells , MPI_LONG_LONG_INT , MPI_SUM , 0 , MPI_COMM_WORLD ) != MPI_SUCCESS )
#endif
       error( "Failed to allreduce edge weights." );
443
444
445
446
        
    /* As of here, only one node needs to compute the partition. */
    if ( nodeID == 0 ) {
    
447
        /* Check that the edge weights are fully symmetric. */
448
        /* for ( cid = 0 ; cid < nr_cells ; cid++ )
449
450
451
452
453
            for ( k = 0 ; k < 26 ; k++ ) {
                cjd = inds[ cid*26 + k ];
                for ( j = 26*cjd ; inds[j] != cid ; j++ );
                if ( weights_e[ cid*26+k ] != weights_e[ j ] )
                    error( "Unsymmetric edge weights detected (%i vs %i)." , weights_e[ cid*26+k ] , weights_e[ j ] );
454
                } */
455
456
457
        /* int w_min = weights_e[0], w_max = weights_e[0], w_tot = weights_e[0];
        for ( k = 1 ; k < 26*nr_cells ; k++ ) {
            w_tot += weights_e[k];
458
459
460
461
            if ( weights_e[k] < w_min )
                w_min = weights_e[k];
            else if ( weights_e[k] > w_max )
                w_max = weights_e[k];
462
463
464
465
466
            }
        message( "edge weights in [ %i , %i ], tot=%i." , w_min , w_max , w_tot );
        w_min = weights_e[0], w_max = weights_e[0]; w_tot = weights_v[0];
        for ( k = 1 ; k < nr_cells ; k++ ) {
            w_tot += weights_v[k];
467
468
469
470
            if ( weights_v[k] < w_min )
                w_min = weights_v[k];
            else if ( weights_v[k] > w_max )
                w_max = weights_v[k];
471
472
            }
        message( "vertex weights in [ %i , %i ], tot=%i." , w_min , w_max , w_tot ); */
473
474
475
476
                
        /* Make sure there are no zero weights. */
        for ( k = 0 ; k < 26*nr_cells ; k++ )
            if ( weights_e[k] == 0 )
Pedro Gonnet's avatar
typo.    
Pedro Gonnet committed
477
                weights_e[k] = 1;
478
        for ( k = 0 ; k < nr_cells ; k++ )
479
            if ( ( weights_v[k] *= vscale ) == 0 )
480
                weights_v[k] = 1;
481
    
482
483
484
485
486
487
488
489
490
491
492
493
494
495
        /* Allocate and fill the connection array. */
        idx_t *offsets;
        if ( ( offsets = (idx_t *)malloc( sizeof(idx_t) * (nr_cells + 1) ) ) == NULL )
            error( "Failed to allocate offsets buffer." );
        offsets[0] = 0;
        for ( k = 0 ; k < nr_cells ; k++ )
            offsets[k+1] = offsets[k] + 26;
            
        /* Set the METIS options. */
        idx_t options[METIS_NOPTIONS];
        METIS_SetDefaultOptions( options );
        options[ METIS_OPTION_OBJTYPE ] = METIS_OBJTYPE_CUT;
        options[ METIS_OPTION_NUMBERING ] = 0;
        options[ METIS_OPTION_CONTIG ] = 1;
496
497
498
        options[ METIS_OPTION_NCUTS ] = 10;
        options[ METIS_OPTION_NITER ] = 20;
        // options[ METIS_OPTION_UFACTOR ] = 1;
499
        
500
501
502
503
        /* Set the initial partition, although this is probably ignored. */
        for ( k = 0 ; k < nr_cells ; k++ )
            nodeIDs[k] = cells[k].nodeID;
            
504
        /* Call METIS. */
505
        idx_t one = 1, idx_nr_cells = nr_cells, idx_nr_nodes = nr_nodes;
506
        idx_t objval;
507
        if ( METIS_PartGraphRecursive( &idx_nr_cells , &one , offsets , inds , weights_v , NULL , weights_e , &idx_nr_nodes , NULL , NULL , options , &objval , nodeIDs ) != METIS_OK )
508
            error( "Call to METIS_PartGraphKway failed." );
509

510
511
512
        /* Dump the 3d array of cell IDs. */
        printf( "engine_repartition: nodeIDs = reshape( [" );
        for ( i = 0 ; i < cdim[0]*cdim[1]*cdim[2] ; i++ )
Pedro Gonnet's avatar
Pedro Gonnet committed
513
            printf( "%i " , (int)nodeIDs[ i ] );
514
        printf("] ,%i,%i,%i);\n",cdim[0],cdim[1],cdim[2]);
515
516
517
518
519
520
521
522
    
        }
        
    /* Broadcast the result of the partition. */
    if ( MPI_Bcast( nodeIDs , nr_cells , MPI_INT , 0 , MPI_COMM_WORLD ) != MPI_SUCCESS )
        error( "Failed to bcast the node IDs." );
        
    /* Set the cell nodeIDs and clear any non-local parts. */
523
    for ( k = 0 ; k < nr_cells ; k++ ) {
524
        cells[k].nodeID = nodeIDs[k];
525
526
527
        if ( nodeIDs[k] == nodeID )
            my_cells += 1;
        }
528
529
530
531
532
533
534
535
536
537
538
539
540
        
    /* Clean up. */
    free( inds );
    free( weights_v );
    free( weights_e );
    free( nodeIDs );
        
    /* Now comes the tricky part: Exchange particles between all nodes.
       This is done in two steps, first allreducing a matrix of 
       how many particles go from where to where, then re-allocating
       the parts array, and emiting the sends and receives.
       Finally, the space, tasks, and proxies need to be rebuilt. */
       
541
542
    /* Redistribute the particles between the nodes. */
    engine_redistribute( e );
543
544
545
546
547
548
549
550
551
552
553
554
        
    /* Make the proxies. */
    engine_makeproxies( e );
        
    /* Tell the engine it should re-build whenever possible */
    e->forcerebuild = 1;
    
#else
    error( "SWIFT was not compiled with MPI and METIS support." );
#endif

    }
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
    
    
/**
 * @brief Add up/down gravity tasks to a cell hierarchy.
 *
 * @param e The #engine.
 * @param c The #cell
 * @param up The upward gravity #task.
 * @param down The downward gravity #task.
 */
 
void engine_addtasks_grav ( struct engine *e , struct cell *c , struct task *up , struct task *down ) {

    /* Link the tasks to this cell. */
    c->grav_up = up;
    c->grav_down = down;
    
    /* Recurse? */
    if ( c->split )
        for ( int k = 0 ; k < 8 ; k++ )
            if ( c->progeny[k] != NULL )
                engine_addtasks_grav( e , c->progeny[k] , up , down );

    }
579
580


581
582
583
584
/**
 * @brief Add send tasks to a hierarchy of cells.
 *
 * @param e The #engine.
585
586
 * @param ci The sending #cell.
 * @param cj The receiving #cell
587
588
589
590
 */

void engine_addtasks_send ( struct engine *e , struct cell *ci , struct cell *cj ) {

Pedro Gonnet's avatar
Pedro Gonnet committed
591
    int k;
592
    struct link *l = NULL;
593
    struct scheduler *s = &e->sched;
594
595

    /* Check if any of the density tasks are for the target node. */
596
597
598
    for ( l = ci->density ; l != NULL ; l = l->next )
        if ( l->t->ci->nodeID == cj->nodeID ||
             ( l->t->cj != NULL && l->t->cj->nodeID == cj->nodeID ) )
599
600
601
            break;

    /* If so, attach send tasks. */
602
    if ( l != NULL ) {
603
604

        /* Create the tasks. */
605
606
        struct task *t_xv = scheduler_addtask( &e->sched , task_type_send , task_subtype_none , 2*ci->tag , 0 , ci , cj , 0 );
        struct task *t_rho = scheduler_addtask( &e->sched , task_type_send , task_subtype_none , 2*ci->tag + 1 , 0 , ci , cj , 0 );
607
608

        /* The send_rho task depends on the cell's ghost task. */
609
        scheduler_addunlock( s , ci->super->ghost , t_rho );
610
611

        /* The send_rho task should unlock the super-cell's kick2 task. */
612
        scheduler_addunlock( s , t_rho , ci->super->kick2 );
613
614

        /* The send_xv task should unlock the super-cell's ghost task. */
615
        scheduler_addunlock( s , t_xv , ci->super->ghost );
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638

        }
        
    /* Recurse? */
    else if ( ci->split )
        for ( k = 0 ; k < 8 ; k++ )
            if ( ci->progeny[k] != NULL )
                engine_addtasks_send( e , ci->progeny[k] , cj );

    }


/**
 * @brief Add recv tasks to a hierarchy of cells.
 *
 * @param e The #engine.
 * @param c The #cell.
 * @param t_xv The recv_xv #task, if it has already been created.
 * @param t_rho The recv_rho #task, if it has already been created.
 */

void engine_addtasks_recv ( struct engine *e , struct cell *c , struct task *t_xv , struct task *t_rho ) {

Pedro Gonnet's avatar
Pedro Gonnet committed
639
    int k;
640
    struct scheduler *s = &e->sched;
641
642

    /* Do we need to construct a recv task? */
643
    if ( t_xv == NULL && c->nr_density > 0 ) {
644
645
    
        /* Create the tasks. */
646
647
        t_xv = c->recv_xv = scheduler_addtask( &e->sched , task_type_recv , task_subtype_none , 2*c->tag , 0 , c , NULL , 0 );
        t_rho = c->recv_rho = scheduler_addtask( &e->sched , task_type_recv , task_subtype_none , 2*c->tag + 1 , 0 , c , NULL , 0 );
648
649
650
        
        }
        
651
652
653
654
655
656
657
658
659
660
    /* Add dependencies. */
    for ( struct link *l = c->density ; l != NULL ; l = l->next ) {
        scheduler_addunlock( s , t_xv , l->t );
        scheduler_addunlock( s , l->t , t_rho );
        }
    for ( struct link *l = c->force ; l != NULL ; l = l->next )
        scheduler_addunlock( s , t_rho , l->t );
    if ( c->sorts != NULL )
        scheduler_addunlock( s , t_xv , c->sorts );
    
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
    /* Recurse? */
    if ( c->split )
        for ( k = 0 ; k < 8 ; k++ )
            if ( c->progeny[k] != NULL )
                engine_addtasks_recv( e , c->progeny[k] , t_xv , t_rho );

    }


/**
 * @brief Exchange cell structures with other nodes.
 *
 * @param e The #engine.
 */
 
void engine_exchange_cells ( struct engine *e ) {

#ifdef WITH_MPI

    int j, k, pid, count = 0;
    struct pcell *pcells;
682
683
684
685
    struct space *s = e->s;
    struct cell *cells = s->cells;
    int nr_cells = s->nr_cells;
    int nr_proxies = e->nr_proxies;
686
    int offset[ nr_cells ];
687
688
    MPI_Request reqs_in[ engine_maxproxies ];
    MPI_Request reqs_out[ engine_maxproxies ];
689
    MPI_Status status;
690
    struct part *parts = &s->parts[ s->nr_parts ];
691
692
693
694
695
696
697
698
699
700
701
702
703
    
    /* Run through the cells and get the size of the ones that will be sent off. */
    for ( k = 0 ; k < nr_cells ; k++ ) {
        offset[k] = count;
        if ( cells[k].sendto )
            count += ( cells[k].pcell_size = cell_getsize( &cells[k] ) );
        }
        
    /* Allocate the pcells. */
    if ( ( pcells = (struct pcell *)malloc( sizeof(struct pcell) * count ) ) == NULL )
        error( "Failed to allocate pcell buffer." );
        
    /* Pack the cells. */
704
    cell_next_tag = 0;
705
706
707
708
709
710
711
    for ( k = 0 ; k < nr_cells ; k++ )
        if ( cells[k].sendto ) {
            cell_pack( &cells[k] , &pcells[ offset[k] ] );
            cells[k].pcell = &pcells[ offset[k] ];
            }

    /* Launch the proxies. */
712
    for ( k = 0 ; k < nr_proxies ; k++ ) {
713
        proxy_cells_exch1( &e->proxies[k] );
714
        reqs_in[k] = e->proxies[k].req_cells_count_in;
Pedro Gonnet's avatar
Pedro Gonnet committed
715
        reqs_out[k] = e->proxies[k].req_cells_count_out;
716
717
718
        }
        
    /* Wait for each count to come in and start the recv. */
719
720
    for ( k = 0 ; k < nr_proxies ; k++ ) {
        if ( MPI_Waitany( nr_proxies , reqs_in , &pid , &status ) != MPI_SUCCESS ||
721
722
723
724
725
726
             pid == MPI_UNDEFINED )
            error( "MPI_Waitany failed." );
        // message( "request from proxy %i has arrived." , pid );
        proxy_cells_exch2( &e->proxies[pid] );
        }
        
Pedro Gonnet's avatar
Pedro Gonnet committed
727
728
729
730
    /* Wait for all the sends to have finnished too. */
    if ( MPI_Waitall( nr_proxies , reqs_out , MPI_STATUSES_IGNORE ) != MPI_SUCCESS )
        error( "MPI_Waitall on sends failed." );
        
731
    /* Set the requests for the cells. */
732
733
734
735
    for ( k = 0 ; k < nr_proxies ; k++ ) {
        reqs_in[k] = e->proxies[k].req_cells_in;
        reqs_out[k] = e->proxies[k].req_cells_out;
        }
736
737
    
    /* Wait for each pcell array to come in from the proxies. */
738
739
    for ( k = 0 ; k < nr_proxies ; k++ ) {
        if ( MPI_Waitany( nr_proxies , reqs_in , &pid , &status ) != MPI_SUCCESS ||
740
741
             pid == MPI_UNDEFINED )
            error( "MPI_Waitany failed." );
742
743
744
745
746
747
        // message( "cell data from proxy %i has arrived." , pid );
        for ( count = 0 , j = 0 ; j < e->proxies[pid].nr_cells_in ; j++ )
            count += cell_unpack( &e->proxies[pid].pcells_in[count] , e->proxies[pid].cells_in[j] , e->s );
        }
        
    /* Wait for all the sends to have finnished too. */
748
    if ( MPI_Waitall( nr_proxies , reqs_out , MPI_STATUSES_IGNORE ) != MPI_SUCCESS )
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
        error( "MPI_Waitall on sends failed." );
        
    /* Count the number of particles we need to import and re-allocate
       the buffer if needed. */
    for ( count = 0 , k = 0 ; k < nr_proxies ; k++ )
        for ( j = 0 ; j < e->proxies[k].nr_cells_in ; j++ )
            count += e->proxies[k].cells_in[j]->count;
    if ( count > s->size_parts_foreign ) {
        if ( s->parts_foreign != NULL )
            free( s->parts_foreign );
        s->size_parts_foreign = 1.1 * count;
        if ( posix_memalign( (void **)&s->parts_foreign , part_align , sizeof(struct part) * s->size_parts_foreign ) != 0 )
            error( "Failed to allocate foreign part data." );
        }
        
    /* Unpack the cells and link to the particle data. */
    parts = s->parts_foreign;
    for ( k = 0 ; k < nr_proxies ; k++ ) {
        for ( count = 0 , j = 0 ; j < e->proxies[k].nr_cells_in ; j++ ) {
            count += cell_link( e->proxies[k].cells_in[j] , parts );
            parts = &parts[ e->proxies[k].cells_in[j]->count ];
770
771
            }
        }
772
773
774
775
776
    s->nr_parts_foreign = parts - s->parts_foreign;
        
    /* Is the parts buffer large enough? */
    if ( s->nr_parts_foreign > s->size_parts_foreign )
        error( "Foreign parts buffer too small." );
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
        
    /* Free the pcell buffer. */
    free( pcells );
    
#else
    error( "SWIFT was not compiled with MPI support." );
#endif

    }


/**
 * @brief Exchange straying parts with other nodes.
 *
 * @param e The #engine.
 * @param parts An array of straying parts.
 * @param xparts The corresponding xparts.
 * @param ind The ID of the foreign #cell.
 * @param N The number of stray parts.
 *
 * @return The number of arrived parts copied to parts and xparts.
 */
 
int engine_exchange_strays ( struct engine *e , struct part *parts , struct xpart *xparts , int *ind , int N ) {

#ifdef WITH_MPI

804
    int k, pid, count = 0, nr_in = 0, nr_out = 0;
Pedro Gonnet's avatar
Pedro Gonnet committed
805
806
    MPI_Request reqs_in[ 2*engine_maxproxies ];
    MPI_Request reqs_out[ 2*engine_maxproxies ];
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
    MPI_Status status;
    struct proxy *p;

    /* Re-set the proxies. */
    for ( k = 0 ; k < e->nr_proxies ; k++ )
        e->proxies[k].nr_parts_out = 0;
    
    /* Put the parts into the corresponding proxies. */
    for ( k = 0 ; k < N ; k++ ) {
        pid = e->proxy_ind[ e->s->cells[ ind[k] ].nodeID ];
        if ( pid < 0 )
            error( "Do not have a proxy for the requested nodeID." );
        proxy_parts_load( &e->proxies[pid] , &parts[k] , &xparts[k] , 1 );
        }
    
    /* Launch the proxies. */
    for ( k = 0 ; k < e->nr_proxies ; k++ ) {
        proxy_parts_exch1( &e->proxies[k] );
825
        reqs_in[k] = e->proxies[k].req_parts_count_in;
Pedro Gonnet's avatar
Pedro Gonnet committed
826
        reqs_out[k] = e->proxies[k].req_parts_count_out;
827
828
829
830
        }
        
    /* Wait for each count to come in and start the recv. */
    for ( k = 0 ; k < e->nr_proxies ; k++ ) {
831
        if ( MPI_Waitany( e->nr_proxies , reqs_in , &pid , &status ) != MPI_SUCCESS ||
832
833
834
835
836
837
             pid == MPI_UNDEFINED )
            error( "MPI_Waitany failed." );
        // message( "request from proxy %i has arrived." , pid );
        proxy_parts_exch2( &e->proxies[pid] );
        }
        
Pedro Gonnet's avatar
Pedro Gonnet committed
838
839
840
841
842
    /* Wait for all the sends to have finnished too. */
    if ( MPI_Waitall( e->nr_proxies , reqs_out , MPI_STATUSES_IGNORE ) != MPI_SUCCESS )
        error( "MPI_Waitall on sends failed." );
        
    /* Return the number of harvested parts. */
843
    /* Set the requests for the particle data. */
844
    for ( k = 0 ; k < e->nr_proxies ; k++ ) {
845
        if ( e->proxies[k].nr_parts_in > 0 ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
846
847
            reqs_in[2*k] = e->proxies[k].req_parts_in;
            reqs_in[2*k+1] = e->proxies[k].req_xparts_in;
848
849
850
            nr_in += 1;
            }
        else
Pedro Gonnet's avatar
Pedro Gonnet committed
851
            reqs_in[2*k] = reqs_in[2*k+1] = MPI_REQUEST_NULL;
852
        if ( e->proxies[k].nr_parts_out > 0 ) {
Pedro Gonnet's avatar
Pedro Gonnet committed
853
854
            reqs_out[2*k] = e->proxies[k].req_parts_out;
            reqs_out[2*k+1] = e->proxies[k].req_xparts_out;
855
856
857
            nr_out += 1;
            }
        else
Pedro Gonnet's avatar
Pedro Gonnet committed
858
            reqs_out[2*k] = reqs_out[2*k+1] = MPI_REQUEST_NULL;
859
        }
860
861
862
    
    /* Wait for each part array to come in and collect the new
       parts from the proxies. */
863
    for ( k = 0 ; k < 2*(nr_in + nr_out) ; k++ ) {
864
865
        int err;
        if ( ( err = MPI_Waitany( 2*e->nr_proxies , reqs_in , &pid , &status ) ) != MPI_SUCCESS ) {
866
867
868
        char buff[ MPI_MAX_ERROR_STRING ];
        int res;
        MPI_Error_string( err , buff , &res );
869
            error( "MPI_Waitany failed (%s)." , buff );
870
871
872
        }
    if ( pid == MPI_UNDEFINED )
        break;
873
874
        if ( pid == MPI_UNDEFINED )
            break;
875
        // message( "request from proxy %i has arrived." , pid );
Pedro Gonnet's avatar
Pedro Gonnet committed
876
877
878
879
880
881
882
883
884
885
886
        if ( reqs_in[pid & ~1] == MPI_REQUEST_NULL &&
             reqs_in[pid | 1 ] == MPI_REQUEST_NULL ) {
            p = &e->proxies[pid/2];
            memcpy( &parts[count] , p->parts_in , sizeof(struct part) * p->nr_parts_in );
            memcpy( &xparts[count] , p->xparts_in , sizeof(struct xpart) * p->nr_parts_in );
            count += p->nr_parts_in;
            /* for ( int k = 0 ; k < p->nr_parts_in ; k++ )
                message( "received particle %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i." ,
                    p->parts_in[k].id , p->parts_in[k].x[0] , p->parts_in[k].x[1] , p->parts_in[k].x[2] ,
                    p->parts_in[k].h , p->nodeID ); */
            }
887
888
        }
    
889
    /* Wait for all the sends to have finnished too. */
890
    if ( nr_out > 0 )
Pedro Gonnet's avatar
Pedro Gonnet committed
891
        if ( MPI_Waitall( 2*e->nr_proxies , reqs_out , MPI_STATUSES_IGNORE ) != MPI_SUCCESS )
892
            error( "MPI_Waitall on sends failed." );
893
        
894
895
896
897
898
899
900
901
902
903
904
    /* Return the number of harvested parts. */
    return count;
    
#else
    error( "SWIFT was not compiled with MPI support." );
    return 0;
#endif

    }


905
/**
906
 * @brief Fill the #space's task list.
907
 *
908
 * @param e The #engine we are working with.
909
910
 */
 
911
void engine_maketasks ( struct engine *e ) {
912
913

    struct space *s = e->s;
914
    struct scheduler *sched = &e->sched;
915
916
917
    struct cell *cells = s->cells;
    int nr_cells = s->nr_cells;
    int nodeID = e->nodeID;
918
919
920
921
922
923
    int i, j, k, ii, jj, kk, iii, jjj, kkk, cid, cjd, sid;
    int *cdim = s->cdim;
    struct task *t, *t2;
    struct cell *ci, *cj;

    /* Re-set the scheduler. */
924
    scheduler_reset( sched , s->tot_cells * engine_maxtaskspercell );
925
926
927
928
929
930
    
    /* Run through the highest level of cells and add pairs. */
    for ( i = 0 ; i < cdim[0] ; i++ )
        for ( j = 0 ; j < cdim[1] ; j++ )
            for ( k = 0 ; k < cdim[2] ; k++ ) {
                cid = cell_getid( cdim , i , j , k );
931
                if ( cells[cid].count == 0 )
932
                    continue;
933
                ci = &cells[cid];
934
935
                if ( ci->count == 0 )
                    continue;
936
                if ( ci->nodeID == nodeID )
937
                    scheduler_addtask( sched , task_type_self , task_subtype_density , 0 , 0 , ci , NULL , 0 );
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
                for ( ii = -1 ; ii < 2 ; ii++ ) {
                    iii = i + ii;
                    if ( !s->periodic && ( iii < 0 || iii >= cdim[0] ) )
                        continue;
                    iii = ( iii + cdim[0] ) % cdim[0];
                    for ( jj = -1 ; jj < 2 ; jj++ ) {
                        jjj = j + jj;
                        if ( !s->periodic && ( jjj < 0 || jjj >= cdim[1] ) )
                            continue;
                        jjj = ( jjj + cdim[1] ) % cdim[1];
                        for ( kk = -1 ; kk < 2 ; kk++ ) {
                            kkk = k + kk;
                            if ( !s->periodic && ( kkk < 0 || kkk >= cdim[2] ) )
                                continue;
                            kkk = ( kkk + cdim[2] ) % cdim[2];
                            cjd = cell_getid( cdim , iii , jjj , kkk );
954
                            cj = &cells[cjd];
955
                            if ( cid >= cjd || cj->count == 0 || 
956
                                 ( ci->nodeID != nodeID && cj->nodeID != nodeID ) )
957
958
                                continue;
                            sid = sortlistID[ (kk+1) + 3*( (jj+1) + 3*(ii+1) ) ];
959
                            scheduler_addtask( sched , task_type_pair , task_subtype_density , sid , 0 , ci , cj , 1 );
960
961
962
963
                            }
                        }
                    }
                }
964
965
                
    /* Add the gravity mm tasks. */
966
967
968
969
970
971
972
    for ( i = 0 ; i < nr_cells ; i++ )
        if ( cells[i].gcount > 0 ) {
            scheduler_addtask( sched , task_type_grav_mm , task_subtype_none , -1 , 0 , &cells[i] , NULL , 0 );
            for ( j = i+1 ; j < nr_cells ; j++ )
                if ( cells[j].gcount > 0 )
                    scheduler_addtask( sched , task_type_grav_mm , task_subtype_none , -1 , 0 , &cells[i] , &cells[j] , 0 );
            }
973
        
974
975
976
    /* Split the tasks. */
    scheduler_splittasks( sched );
    
977
978
979
    /* Allocate the list of cell-task links. The maximum number of links
       is the number of cells (s->tot_cells) times the number of neighbours (27)
       times the number of interaction types (2, density and force). */
980
981
    if ( e->links != NULL )
        free( e->links );
982
    if ( ( e->links = malloc( sizeof(struct link) * s->tot_cells * 27 * 2 ) ) == NULL )
983
984
985
        error( "Failed to allocate cell-task links." );
    e->nr_links = 0;
    
986
987
    /* Add the gravity up/down tasks at the top-level cells and push them down. */
    for ( k = 0 ; k < nr_cells ; k++ )
988
        if ( cells[k].nodeID == nodeID && cells[k].gcount > 0 ) {
989
990
991
992
993
994
995
996
997
998
        
            /* Create tasks at top level. */
            struct task *up = scheduler_addtask( sched , task_type_grav_up , task_subtype_none , 0 , 0 , &cells[k] , NULL , 0 );
            struct task *down = scheduler_addtask( sched , task_type_grav_down , task_subtype_none , 0 , 0 , &cells[k] , NULL , 0 );
            
            /* Push tasks down the cell hierarchy. */
            engine_addtasks_grav( e , &cells[k] , up , down );
            
            }
    
999
1000
1001
1002
1003
    /* Count the number of tasks associated with each cell and
       store the density tasks in each cell, and make each sort
       depend on the sorts of its progeny. */
    // #pragma omp parallel for private(t,j)
    for ( k = 0 ; k < sched->nr_tasks ; k++ ) {
1004
1005
        
        /* Get the current task. */
1006
1007
1008
        t = &sched->tasks[k];
        if ( t->skip )
            continue;
1009
1010
            
        /* Link sort tasks together. */
1011
        if ( t->type == task_type_sort && t->ci->split )
1012
1013
            for ( j = 0 ; j < 8 ; j++ )
                if ( t->ci->progeny[j] != NULL && t->ci->progeny[j]->sorts != NULL ) {
1014
                    t->ci->progeny[j]->sorts->skip = 0;
1015
                    scheduler_addunlock( sched , t->ci->progeny[j]->sorts , t );
1016
                    }
1017
1018
                    
        /* Link density tasks to cells. */
1019
1020
1021
        if ( t->type == task_type_self ) {
            atomic_inc( &t->ci->nr_tasks );
            if ( t->subtype == task_subtype_density ) {
1022
1023
                t->ci->density = engine_addlink( e , t->ci->density , t );
                atomic_inc( &t->ci->nr_density );
1024
1025
1026
1027
1028
1029
                }
            }
        else if ( t->type == task_type_pair ) {
            atomic_inc( &t->ci->nr_tasks );
            atomic_inc( &t->cj->nr_tasks );
            if ( t->subtype == task_subtype_density ) {
1030
1031
1032
1033
                t->ci->density = engine_addlink( e , t->ci->density , t );
                atomic_inc( &t->ci->nr_density );
                t->cj->density = engine_addlink( e , t->cj->density , t );
                atomic_inc( &t->cj->nr_density );
1034
1035
1036
1037
1038
1039
1040
                }
            }
        else if ( t->type == task_type_sub ) {
            atomic_inc( &t->ci->nr_tasks );
            if ( t->cj != NULL )
                atomic_inc( &t->cj->nr_tasks );
            if ( t->subtype == task_subtype_density ) {
1041
1042
1043
1044
1045
                t->ci->density = engine_addlink( e , t->ci->density , t );
                atomic_inc( &t->ci->nr_density );
                if ( t->cj != NULL ) {
                    t->cj->density = engine_addlink( e , t->cj->density , t );
                    atomic_inc( &t->cj->nr_density );
1046
1047
                    }
                }
1048
            }
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
            
        /* Link gravity multipole tasks to the up/down tasks. */
        if ( t->type == task_type_grav_mm ||
             ( t->type == task_type_sub && t->subtype == task_subtype_grav ) ) {
            atomic_inc( &t->ci->nr_tasks );
            scheduler_addunlock( sched , t->ci->grav_up , t );
            scheduler_addunlock( sched , t , t->ci->grav_down );
            if ( t->cj != NULL && t->ci->grav_up != t->cj->grav_up ) {
                scheduler_addunlock( sched , t->cj->grav_up , t );
                scheduler_addunlock( sched , t , t->cj->grav_down );
1059
1060
                }
            }
1061
            
1062
1063
        }
        
1064
1065
1066
1067
1068
    /* Append a ghost task to each cell, and add kick2 tasks to the
       super cells. */
    for ( k = 0 ; k < nr_cells ; k++ )
        engine_mkghosts( e , &cells[k] , NULL );
        
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
    /* Run through the tasks and make force tasks for each density task.
       Each force task depends on the cell ghosts and unlocks the kick2 task
       of its super-cell. */
    kk = sched->nr_tasks;
    // #pragma omp parallel for private(t,t2)
    for ( k = 0 ; k < kk ; k++ ) {
    
        /* Get a pointer to the task. */
        t = &sched->tasks[k];
        
        /* Skip? */
        if ( t->skip )
            continue;
1082
            
1083
1084
        /* Self-interaction? */
        if ( t->type == task_type_self && t->subtype == task_subtype_density ) {
1085
            scheduler_addunlock( sched , t , t->ci->super->ghost );
1086
            t2 = scheduler_addtask( sched , task_type_self , task_subtype_force , 0 , 0 , t->ci , NULL , 0 );
1087
            scheduler_addunlock( sched , t->ci->super->ghost , t2 );
1088
            scheduler_addunlock( sched , t2 , t->ci->super->kick2 );
1089
1090
            t->ci->force = engine_addlink( e , t->ci->force , t2 );
            atomic_inc( &t->ci->nr_force );
1091
1092
1093
1094
1095
            }
            
        /* Otherwise, pair interaction? */
        else if ( t->type == task_type_pair && t->subtype == task_subtype_density ) {
            t2 = scheduler_addtask( sched , task_type_pair , task_subtype_force , 0 , 0 , t->ci , t->cj , 0 );
1096
            if ( t->ci->nodeID == nodeID ) {
1097
                scheduler_addunlock( sched , t , t->ci->super->ghost );
1098
                scheduler_addunlock( sched , t->ci->super->ghost , t2 );
1099
                scheduler_addunlock( sched , t2 , t->ci->super->kick2 );
1100
                }
1101
            if ( t->cj->nodeID == nodeID && t->ci->super != t->cj->super ) {
1102
1103
1104
                scheduler_addunlock( sched , t , t->cj->super->ghost );
                scheduler_addunlock( sched , t->cj->super->ghost , t2 );
                scheduler_addunlock( sched , t2 , t->cj->super->kick2 );
1105
                }
1106
1107
1108
1109
            t->ci->force = engine_addlink( e , t->ci->force , t2 );
            atomic_inc( &t->ci->nr_force );
            t->cj->force = engine_addlink( e , t->cj->force , t2 );
            atomic_inc( &t->cj->nr_force );
1110
1111
1112
1113
1114
            }
    
        /* Otherwise, sub interaction? */
        else if ( t->type == task_type_sub && t->subtype == task_subtype_density ) {
            t2 = scheduler_addtask( sched , task_type_sub , task_subtype_force , t->flags , 0 , t->ci , t->cj , 0 );
1115
            if ( t->ci->nodeID == nodeID ) {
1116
                scheduler_addunlock( sched , t , t->ci->super->ghost );
1117
                scheduler_addunlock( sched , t->ci->super->ghost , t2 );
1118
                scheduler_addunlock( sched , t2 , t->ci->super->kick2 );
1119
                }
1120
            if ( t->cj != NULL && t->cj->nodeID == nodeID && t->ci->super != t->cj->super ) {
1121
1122
1123
1124
1125
1126
1127
1128
1129
                scheduler_addunlock( sched , t , t->cj->super->ghost );
                scheduler_addunlock( sched , t->cj->super->ghost , t2 );
                scheduler_addunlock( sched , t2 , t->cj->super->kick2 );
                }
            t->ci->force = engine_addlink( e , t->ci->force , t2 );
            atomic_inc( &t->ci->nr_force );
            if ( t->cj != NULL ) {
                t->cj->force = engine_addlink( e , t->cj->force , t2 );
                atomic_inc( &t->cj->nr_force );
1130
                }
1131
1132
            }
            
1133
        /* Kick2 tasks should rely on the grav_down tasks of their cell. */
1134
        else if ( t->type == task_type_kick2 && t->ci->grav_down != NULL )
1135
1136
            scheduler_addunlock( sched , t->ci->grav_down , t );
            
1137
1138
        }
        
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
    /* Add the communication tasks if MPI is being used. */
    #ifdef WITH_MPI
        
        /* Loop over the proxies. */
        for ( int pid = 0 ; pid < e->nr_proxies ; pid++ ) {
        
            /* Get a handle on the proxy. */
            struct proxy *p = &e->proxies[pid];
            
            /* Loop through the proxy's incomming cells and add the
               recv tasks. */
            for ( k = 0 ; k < p->nr_cells_in ; k++ )
                engine_addtasks_recv( e , p->cells_in[k] , NULL , NULL );
            
            /* Loop through the proxy's outgoing cells and add the
               send tasks. */
1155
            for ( k = 0 ; k < p->nr_cells_out ; k++ )
1156
1157
1158
1159
1160
1161
                engine_addtasks_send( e , p->cells_out[k] , p->cells_in[0] );
            
            }
        
    #endif
        
1162
1163
    /* Rank the tasks. */
    scheduler_ranktasks( sched );
Pedro Gonnet's avatar
Pedro Gonnet committed
1164
    
1165
1166
1167
    /* Weight the tasks. */
    scheduler_reweight( sched );
            
Pedro Gonnet's avatar
Pedro Gonnet committed
1168
1169
1170
    /* Set the tasks age. */
    e->tasks_age = 0;
            
1171
1172
    }
    
1173
    
1174

1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
/**
 * @brief Mark tasks to be skipped and set the sort flags accordingly.
 * 
 * @return 1 if the space has to be rebuilt, 0 otherwise.
 */
 
int engine_marktasks ( struct engine *e ) {

    struct scheduler *s = &e->sched;
    int k, nr_tasks = s->nr_tasks, *ind = s->tasks_ind;
    struct task *t, *tasks = s->tasks;
    float dt_step = e->dt_step;
    struct cell *ci, *cj;
1188
    // ticks tic = getticks();
1189
    
1190
1191
    /* Muc less to do here if we're on a fixed time-step. */
    if ( !( e->policy & engine_policy_multistep ) ) {
1192
    
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
        /* Run through the tasks and mark as skip or not. */
        for ( k = 0 ; k < nr_tasks ; k++ ) {

            /* Get a handle on the kth task. */
            t = &tasks[ ind[k] ];

            /* Pair? */
            if ( t->type == task_type_pair || ( t->type == task_type_sub && t->cj != NULL ) ) {

                /* Local pointers. */
                ci = t->ci;
                cj = t->cj;

                /* Too much particle movement? */
                if ( t->tight &&
                     ( fmaxf( ci->h_max , cj->h_max ) + ci->dx_max + cj->dx_max > cj->dmin || 
                       ci->dx_max > space_maxreldx*ci->h_max || cj->dx_max > space_maxreldx*cj->h_max ) )
                    return 1;

                }
1213
1214
1215
1216
1217
1218
1219
1220
1221
                
            /* Sort? */
            else if ( t->type == task_type_sort ) {
            
                /* If all the sorts have been done, make this task implicit. */
                if ( !( t->flags & (t->flags ^ t->ci->sorted ) ) )
                    t->implicit = 1;
            
                }
1222

1223
1224
            }
            
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
        }
    
    else {
    
        /* Run through the tasks and mark as skip or not. */
        for ( k = 0 ; k < nr_tasks ; k++ ) {

            /* Get a handle on the kth task. */
            t = &tasks[ ind[k] ];

            /* Sort-task? Note that due to the task ranking, the sorts
               will all come before the pairs. */
            if ( t->type == task_type_sort ) {

                /* Re-set the flags. */
                t->flags = 0;
                t->skip = 1;

                }

            /* Single-cell task? */
            else if ( t->type == task_type_self ||
                      t->type == task_type_ghost ||
                    ( t->type == task_type_sub && t->cj == NULL ) ) {

                /* Set this task's skip. */
                t->skip = ( t->ci->dt_min > dt_step );

                }

            /* Pair? */
            else if ( t->type == task_type_pair || ( t->type == task_type_sub && t->cj != NULL ) ) {

                /* Local pointers