space.c 45.2 KB
Newer Older
1
2
3
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
*               2016 Peter W. Draper (p.w.draper@durham.ac.uk)
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
*
* 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
20
21
22
23
24
25
26
27

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

/* Some standard headers. */
#include <float.h>
#include <limits.h>
#include <math.h>
28
#include <string.h>
29
#include <stdlib.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
30

31
32
/* MPI headers. */
#ifdef WITH_MPI
33
#include <mpi.h>
34
35
#endif

36
37
38
/* This object's header. */
#include "space.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
39
/* Local headers. */
40
#include "atomic.h"
41
#include "engine.h"
42
#include "error.h"
43
#include "kernel_hydro.h"
44
#include "lock.h"
45
#include "minmax.h"
46
#include "runner.h"
47
#include "tools.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
48

49
50
51
/* Shared sort structure. */
struct parallel_sort space_sort_struct;

Pedro Gonnet's avatar
Pedro Gonnet committed
52
53
/* Split size. */
int space_splitsize = space_splitsize_default;
54
int space_subsize = space_subsize_default;
55
int space_maxsize = space_maxsize_default;
Pedro Gonnet's avatar
Pedro Gonnet committed
56
57
58

/* Map shift vector to sortlist. */
const int sortlistID[27] = {
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
    /* ( -1 , -1 , -1 ) */ 0,
    /* ( -1 , -1 ,  0 ) */ 1,
    /* ( -1 , -1 ,  1 ) */ 2,
    /* ( -1 ,  0 , -1 ) */ 3,
    /* ( -1 ,  0 ,  0 ) */ 4,
    /* ( -1 ,  0 ,  1 ) */ 5,
    /* ( -1 ,  1 , -1 ) */ 6,
    /* ( -1 ,  1 ,  0 ) */ 7,
    /* ( -1 ,  1 ,  1 ) */ 8,
    /* (  0 , -1 , -1 ) */ 9,
    /* (  0 , -1 ,  0 ) */ 10,
    /* (  0 , -1 ,  1 ) */ 11,
    /* (  0 ,  0 , -1 ) */ 12,
    /* (  0 ,  0 ,  0 ) */ 0,
    /* (  0 ,  0 ,  1 ) */ 12,
    /* (  0 ,  1 , -1 ) */ 11,
    /* (  0 ,  1 ,  0 ) */ 10,
    /* (  0 ,  1 ,  1 ) */ 9,
    /* (  1 , -1 , -1 ) */ 8,
    /* (  1 , -1 ,  0 ) */ 7,
    /* (  1 , -1 ,  1 ) */ 6,
    /* (  1 ,  0 , -1 ) */ 5,
    /* (  1 ,  0 ,  0 ) */ 4,
    /* (  1 ,  0 ,  1 ) */ 3,
    /* (  1 ,  1 , -1 ) */ 2,
    /* (  1 ,  1 ,  0 ) */ 1,
    /* (  1 ,  1 ,  1 ) */ 0};

87
88
89
90
91
92
93
94
95
96
97
98
/**
 * @brief Get the shift-id of the given pair of cells, swapping them
 *      if need be.
 *
 * @param s The space
 * @param ci Pointer to first #cell.
 * @param cj Pointer second #cell.
 * @param shift Vector from ci to cj.
 *
 * @return The shift ID and set shift, may or may not swap ci and cj.
 */

99
100
101
102
int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
                 double *shift) {

  /* Get the relative distance between the pairs, wrapping. */
103
104
105
  const int periodic = s->periodic;
  double dx[3];
  for (int k = 0; k < 3; k++) {
106
107
108
109
110
111
112
113
114
115
116
    dx[k] = (*cj)->loc[k] - (*ci)->loc[k];
    if (periodic && dx[k] < -s->dim[k] / 2)
      shift[k] = s->dim[k];
    else if (periodic && dx[k] > s->dim[k] / 2)
      shift[k] = -s->dim[k];
    else
      shift[k] = 0.0;
    dx[k] += shift[k];
  }

  /* Get the sorting index. */
117
  int sid = 0;
118
  for (int k = 0; k < 3; k++)
119
120
121
122
    sid = 3 * sid + ((dx[k] < 0.0) ? 0 : ((dx[k] > 0.0) ? 2 : 1));

  /* Switch the cells around? */
  if (runner_flip[sid]) {
123
    struct cell *temp = *ci;
124
125
    *ci = *cj;
    *cj = temp;
126
    for (int k = 0; k < 3; k++) shift[k] = -shift[k];
127
128
129
130
131
132
  }
  sid = sortlistID[sid];

  /* Return the sort ID. */
  return sid;
}
133

134
/**
135
 * @brief Recursively dismantle a cell tree.
136
137
 *
 */
138
139
140
141

void space_rebuild_recycle(struct space *s, struct cell *c) {

  if (c->split)
142
    for (int k = 0; k < 8; k++)
143
144
145
146
147
148
149
      if (c->progeny[k] != NULL) {
        space_rebuild_recycle(s, c->progeny[k]);
        space_recycle(s, c->progeny[k]);
        c->progeny[k] = NULL;
      }
}

150
/**
151
 * @brief Re-build the cell grid.
152
 *
153
154
 * @param s The #space.
 * @param cell_max Maximum cell edge length.
155
 * @param verbose Print messages to stdout or not.
156
 */
157

158
void space_regrid(struct space *s, double cell_max, int verbose) {
159

160
  const size_t nr_parts = s->nr_parts;
161
  struct cell *restrict c;
162
  ticks tic = getticks();
163
164
165

  /* Run through the parts and get the current h_max. */
  // tic = getticks();
166
  float h_max = s->cell_min / kernel_gamma / space_stretch;
167
  if (s->cells != NULL) {
168
    for (int k = 0; k < s->nr_cells; k++) {
169
      if (s->cells[k].h_max > h_max) h_max = s->cells[k].h_max;
170
    }
171
  } else {
Pedro Gonnet's avatar
Pedro Gonnet committed
172
    for (size_t k = 0; k < nr_parts; k++) {
173
174
175
176
177
178
179
180
181
182
183
184
      if (s->parts[k].h > h_max) h_max = s->parts[k].h;
    }
    s->h_max = h_max;
  }

/* If we are running in parallel, make sure everybody agrees on
   how large the largest cell should be. */
#ifdef WITH_MPI
  {
    float buff;
    if (MPI_Allreduce(&h_max, &buff, 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD) !=
        MPI_SUCCESS)
185
      error("Failed to aggregate the rebuild flag across nodes.");
186
187
188
    h_max = buff;
  }
#endif
189
  if (verbose) message("h_max is %.3e (cell_max=%.3e).", h_max, cell_max);
190
191

  /* Get the new putative cell dimensions. */
192
  int cdim[3];
193
  for (int k = 0; k < 3; k++)
194
195
196
197
198
199
200
201
202
    cdim[k] =
        floor(s->dim[k] / fmax(h_max * kernel_gamma * space_stretch, cell_max));

  /* Check if we have enough cells for periodicity. */
  if (s->periodic && (cdim[0] < 3 || cdim[1] < 3 || cdim[2] < 3))
    error(
        "Must have at least 3 cells in each spatial dimension when periodicity "
        "is switched on.");

203
204
205
  /* In MPI-Land, changing the top-level cell size requires that the
   * global partition is recomputed and the particles redistributed.
   * Be prepared to do that. */
206
#ifdef WITH_MPI
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
  double oldh[3];
  double oldcdim[3];
  int *oldnodeIDs = NULL;
  if (cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] || cdim[2] < s->cdim[2]) {

    /* Capture state of current space. */
    oldcdim[0] = s->cdim[0];
    oldcdim[1] = s->cdim[1];
    oldcdim[2] = s->cdim[2];
    oldh[0] = s->h[0];
    oldh[1] = s->h[1];
    oldh[2] = s->h[2];

    if ((oldnodeIDs = (int *)malloc(sizeof(int) * s->nr_cells)) == NULL)
      error("Failed to allocate temporary nodeIDs.");

    int cid = 0;
    for (int i = 0; i < s->cdim[0]; i++) {
      for (int j = 0; j < s->cdim[1]; j++) {
        for (int k = 0; k < s->cdim[2]; k++) {
          cid = cell_getid(oldcdim, i, j, k);
          oldnodeIDs[cid] = s->cells[cid].nodeID;
        }
      }
    }
  }

234
235
236
237
238
239
240
241
242
#endif

  /* Do we need to re-build the upper-level cells? */
  // tic = getticks();
  if (s->cells == NULL || cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] ||
      cdim[2] < s->cdim[2]) {

    /* Free the old cells, if they were allocated. */
    if (s->cells != NULL) {
243
      for (int k = 0; k < s->nr_cells; k++) {
244
245
246
247
248
249
250
251
        space_rebuild_recycle(s, &s->cells[k]);
        if (s->cells[k].sort != NULL) free(s->cells[k].sort);
      }
      free(s->cells);
      s->maxdepth = 0;
    }

    /* Set the new cell dimensions only if smaller. */
252
    for (int k = 0; k < 3; k++) {
253
254
255
256
      s->cdim[k] = cdim[k];
      s->h[k] = s->dim[k] / cdim[k];
      s->ih[k] = 1.0 / s->h[k];
    }
257
    const float dmin = fminf(s->h[0], fminf(s->h[1], s->h[2]));
258
259
260
261
262
263
264

    /* Allocate the highest level of cells. */
    s->tot_cells = s->nr_cells = cdim[0] * cdim[1] * cdim[2];
    if (posix_memalign((void *)&s->cells, 64,
                       s->nr_cells * sizeof(struct cell)) != 0)
      error("Failed to allocate cells.");
    bzero(s->cells, s->nr_cells * sizeof(struct cell));
265
    for (int k = 0; k < s->nr_cells; k++)
266
267
268
      if (lock_init(&s->cells[k].lock) != 0) error("Failed to init spinlock.");

    /* Set the cell location and sizes. */
269
270
271
    for (int i = 0; i < cdim[0]; i++)
      for (int j = 0; j < cdim[1]; j++)
        for (int k = 0; k < cdim[2]; k++) {
272
273
274
275
276
277
278
279
280
281
282
283
284
          c = &s->cells[cell_getid(cdim, i, j, k)];
          c->loc[0] = i * s->h[0];
          c->loc[1] = j * s->h[1];
          c->loc[2] = k * s->h[2];
          c->h[0] = s->h[0];
          c->h[1] = s->h[1];
          c->h[2] = s->h[2];
          c->dmin = dmin;
          c->depth = 0;
          c->count = 0;
          c->gcount = 0;
          c->super = c;
          lock_init(&c->lock);
Pedro Gonnet's avatar
Pedro Gonnet committed
285
        }
286
287

    /* Be verbose about the change. */
288
289
290
    if (verbose)
      message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1],
              cdim[2]);
291
292
    fflush(stdout);

293
#ifdef WITH_MPI
294
295
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
    if (oldnodeIDs != NULL) {
      /* We have changed the top-level cell dimension, so need to redistribute
       * cells around the nodes. We repartition using the old space node
       * positions as a grid to resample. */
      if (s->e->nodeID == 0)
        message("basic cell dimensions have increased - recalculating the "
                "global partition.");

      if (!partition_space_to_space(oldh, oldcdim, oldnodeIDs, s) ) {

        /* Failed, try another technique that requires no settings. */
        message("Failed to get a new partition, trying less optimal method");
        struct partition initial_partition;
#ifdef HAVE_METIS
        initial_partition.type = INITPART_METIS_NOWEIGHT;
#else
        initial_partition.type = INITPART_VECTORIZE;
#endif
        partition_initial_partition(&initial_partition, s->e->nodeID,
                                    s->e->nr_nodes, s);
      }

      /* Re-distribute the particles to their new nodes. */
      engine_redistribute(s->e);

      /* Make the proxies. */
      engine_makeproxies(s->e);
321

322
323
      /* Finished with these. */
      free(oldnodeIDs);
324
325
    }
#endif
326
  } /* re-build upper-level cells? */
327
328
  // message( "rebuilding upper-level cells took %.3f %s." ,
  // clocks_from_ticks(double)(getticks() - tic), clocks_getunit());
329
330
331
332
333

  /* Otherwise, just clean up the cells. */
  else {

    /* Free the old cells, if they were allocated. */
334
    for (int k = 0; k < s->nr_cells; k++) {
335
336
337
338
339
340
341
342
343
344
345
      space_rebuild_recycle(s, &s->cells[k]);
      s->cells[k].sorts = NULL;
      s->cells[k].nr_tasks = 0;
      s->cells[k].nr_density = 0;
      s->cells[k].nr_force = 0;
      s->cells[k].density = NULL;
      s->cells[k].force = NULL;
      s->cells[k].dx_max = 0.0f;
      s->cells[k].sorted = 0;
      s->cells[k].count = 0;
      s->cells[k].gcount = 0;
Matthieu Schaller's avatar
Matthieu Schaller committed
346
      s->cells[k].init = NULL;
Matthieu Schaller's avatar
Matthieu Schaller committed
347
      s->cells[k].ghost = NULL;
Matthieu Schaller's avatar
Matthieu Schaller committed
348
349
      s->cells[k].drift = NULL;
      s->cells[k].kick = NULL;
350
      s->cells[k].super = &s->cells[k];
351
    }
352
353
    s->maxdepth = 0;
  }
354
355
356
357

  if (verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
358
}
359
360
361
362
363
364

/**
 * @brief Re-build the cells as well as the tasks.
 *
 * @param s The #space in which to update the cells.
 * @param cell_max Maximal cell size.
365
 * @param verbose Print messages to stdout or not
366
367
 *
 */
368

369
void space_rebuild(struct space *s, double cell_max, int verbose) {
370

Matthieu Schaller's avatar
Matthieu Schaller committed
371
  const ticks tic = getticks();
372
373
374
375
376

  /* Be verbose about this. */
  // message( "re)building space..." ); fflush(stdout);

  /* Re-grid if necessary, or just re-set the cell data. */
377
  space_regrid(s, cell_max, verbose);
378

Pedro Gonnet's avatar
Pedro Gonnet committed
379
380
  size_t nr_parts = s->nr_parts;
  size_t nr_gparts = s->nr_gparts;
381
382
  struct cell *restrict cells = s->cells;

Matthieu Schaller's avatar
Matthieu Schaller committed
383
384
385
  const double ih[3] = {s->ih[0], s->ih[1], s->ih[2]};
  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
  const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
386
387
388
389

  /* Run through the particles and get their cell index. */
  // tic = getticks();
  const size_t ind_size = s->size_parts;
390
391
  int *ind;
  if ((ind = (int *)malloc(sizeof(int) * ind_size)) == NULL)
392
    error("Failed to allocate temporary particle indices.");
Pedro Gonnet's avatar
Pedro Gonnet committed
393
  for (size_t k = 0; k < nr_parts; k++) {
394
395
    struct part *restrict p = &s->parts[k];
    for (int j = 0; j < 3; j++)
396
397
398
399
      if (p->x[j] < 0.0)
        p->x[j] += dim[j];
      else if (p->x[j] >= dim[j])
        p->x[j] -= dim[j];
400
    ind[k] =
401
        cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
402
    cells[ind[k]].count++;
403
  }
Pedro Gonnet's avatar
Pedro Gonnet committed
404
405
  // message( "getting particle indices took %.3f %s." ,
  // clocks_from_ticks(getticks() - tic), clocks_getunit()):
406

407
408
409
410
411
412
413
  /* Run through the gravity particles and get their cell index. */
  // tic = getticks();
  const size_t gind_size = s->size_gparts;
  int *gind;
  if ((gind = (int *)malloc(sizeof(int) * gind_size)) == NULL)
    error("Failed to allocate temporary g-particle indices.");
  for (int k = 0; k < nr_gparts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
414
    struct gpart *restrict gp = &s->gparts[k];
415
416
417
418
419
420
421
422
423
424
425
426
    for (int j = 0; j < 3; j++)
      if (gp->x[j] < 0.0)
        gp->x[j] += dim[j];
      else if (gp->x[j] >= dim[j])
        gp->x[j] -= dim[j];
    gind[k] =
        cell_getid(cdim, gp->x[0] * ih[0], gp->x[1] * ih[1], gp->x[2] * ih[2]);
    cells[gind[k]].gcount++;
  }
// message( "getting particle indices took %.3f %s." ,
// clocks_from_ticks(getticks() - tic), clocks_getunit());

427
428
#ifdef WITH_MPI
  /* Move non-local parts to the end of the list. */
429
  const int local_nodeID = s->e->nodeID;
Pedro Gonnet's avatar
Pedro Gonnet committed
430
  for (size_t k = 0; k < nr_parts; k++)
431
    if (cells[ind[k]].nodeID != local_nodeID) {
432
433
      cells[ind[k]].count -= 1;
      nr_parts -= 1;
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
434
      const struct part tp = s->parts[k];
435
436
      s->parts[k] = s->parts[nr_parts];
      s->parts[nr_parts] = tp;
437
438
439
440
441
442
      if (s->parts[k].gpart != NULL) {
        s->parts[k].gpart->part = &s->parts[k];
      }
      if (s->parts[nr_parts].gpart != NULL) {
        s->parts[nr_parts].gpart->part = &s->parts[nr_parts];
      }
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
443
      const struct xpart txp = s->xparts[k];
444
445
      s->xparts[k] = s->xparts[nr_parts];
      s->xparts[nr_parts] = txp;
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
446
      const int t = ind[k];
447
448
      ind[k] = ind[nr_parts];
      ind[nr_parts] = t;
449
450
    }

451
452
  /* Move non-local gparts to the end of the list. */
  for (int k = 0; k < nr_gparts; k++)
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
453
454
    if (cells[gind[k]].nodeID != local_nodeID) {
      cells[gind[k]].gcount -= 1;
455
      nr_gparts -= 1;
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
456
      const struct gpart tp = s->gparts[k];
457
458
      s->gparts[k] = s->gparts[nr_gparts];
      s->gparts[nr_gparts] = tp;
459
460
461
462
463
464
      if (s->gparts[k].id > 0) {
        s->gparts[k].part->gpart = &s->gparts[k];
      }
      if (s->gparts[nr_gparts].id > 0) {
        s->gparts[nr_gparts].part->gpart = &s->gparts[nr_gparts];
      }
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
465
466
467
      const int t = gind[k];
      gind[k] = gind[nr_gparts];
      gind[nr_gparts] = t;
468
469
    }

470
471
  /* Exchange the strays, note that this potentially re-allocates
     the parts arrays. */
472
473
474
  /* TODO: This function also exchanges gparts, but this is shorted-out
     until they are fully implemented. */
  size_t nr_parts_exchanged = s->nr_parts - nr_parts;
475
  size_t nr_gparts_exchanged = s->nr_gparts - nr_gparts;
Pedro Gonnet's avatar
Pedro Gonnet committed
476
477
478
  engine_exchange_strays(s->e, nr_parts, &ind[nr_parts], &nr_parts_exchanged,
                         nr_gparts, &gind[nr_gparts], &nr_gparts_exchanged);

479
  /* Add post-processing, i.e. re-linking/creating of gparts here. */
Pedro Gonnet's avatar
Pedro Gonnet committed
480
481

  /* Set the new particle counts. */
482
  s->nr_parts = nr_parts + nr_parts_exchanged;
483
  s->nr_gparts = nr_gparts + nr_gparts_exchanged;
484
485

  /* Re-allocate the index array if needed.. */
486
  if (s->nr_parts > ind_size) {
487
488
    int *ind_new;
    if ((ind_new = (int *)malloc(sizeof(int) * s->nr_parts)) == NULL)
489
      error("Failed to allocate temporary particle indices.");
490
    memcpy(ind_new, ind, sizeof(size_t) * nr_parts);
491
492
    free(ind);
    ind = ind_new;
493
494
495
  }

  /* Assign each particle to its cell. */
Pedro Gonnet's avatar
Pedro Gonnet committed
496
  for (size_t k = nr_parts; k < s->nr_parts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
497
    const struct part *const p = &s->parts[k];
498
    ind[k] =
499
        cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
500
501
502
503
    cells[ind[k]].count += 1;
    /* if ( cells[ ind[k] ].nodeID != nodeID )
        error( "Received part that does not belong to me (nodeID=%i)." , cells[
       ind[k] ].nodeID ); */
504
  }
505
  nr_parts = s->nr_parts;
506
507
508
#endif

  /* Sort the parts according to their cells. */
509
  space_parts_sort(s, ind, nr_parts, 0, s->nr_cells - 1, verbose);
510
511

  /* Re-link the gparts. */
Pedro Gonnet's avatar
Pedro Gonnet committed
512
  for (size_t k = 0; k < nr_parts; k++)
513
    if (s->parts[k].gpart != NULL) s->parts[k].gpart->part = &s->parts[k];
514

515
  /* Verify space_sort_struct. */
516
  /* for ( k = 1 ; k < nr_parts ; k++ ) {
517
      if ( ind[k-1] > ind[k] ) {
518
519
          error( "Sort failed!" );
          }
520
      else if ( ind[k] != cell_getid( cdim , parts[k].x[0]*ih[0] ,
521
522
523
524
525
     parts[k].x[1]*ih[1] , parts[k].x[2]*ih[2] ) )
          error( "Incorrect indices!" );
      } */

  /* We no longer need the indices as of here. */
526
  free(ind);
527

528
529
530
531
#ifdef WITH_MPI

  /* Re-allocate the index array if needed.. */
  if (s->nr_gparts > gind_size) {
532
533
    int *gind_new;
    if ((gind_new = (int *)malloc(sizeof(int) * s->nr_gparts)) == NULL)
534
      error("Failed to allocate temporary g-particle indices.");
535
    memcpy(gind_new, gind, sizeof(int) * nr_gparts);
536
537
538
539
540
    free(gind);
    gind = gind_new;
  }

  /* Assign each particle to its cell. */
541
  for (int k = nr_gparts; k < s->nr_gparts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
542
    const struct gpart *const p = &s->gparts[k];
543
544
    gind[k] =
        cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
Matthieu Schaller's avatar
Typo    
Matthieu Schaller committed
545
    cells[gind[k]].gcount += 1;
546
547
548
549
550
    /* if ( cells[ ind[k] ].nodeID != nodeID )
        error( "Received part that does not belong to me (nodeID=%i)." , cells[
       ind[k] ].nodeID ); */
  }
  nr_gparts = s->nr_gparts;
551

552
#endif
553
554

  /* Sort the parts according to their cells. */
Matthieu Schaller's avatar
Matthieu Schaller committed
555
  space_gparts_sort(s, gind, nr_gparts, 0, s->nr_cells - 1, verbose);
556
557

  /* Re-link the parts. */
558
  for (int k = 0; k < nr_gparts; k++)
559
    if (s->gparts[k].id > 0) s->gparts[k].part->gpart = &s->gparts[k];
560
561

  /* We no longer need the indices as of here. */
562
  free(gind);
563

564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
  /* Verify that the links are correct */
  /* MATTHIEU: To be commented out once we are happy */
  for (size_t k = 0; k < nr_gparts; ++k) {

    if (s->gparts[k].id > 0) {

      if (s->gparts[k].part->gpart != &s->gparts[k]) error("Linking problem !");

      if (s->gparts[k].x[0] != s->gparts[k].part->x[0] ||
          s->gparts[k].x[1] != s->gparts[k].part->x[1] ||
          s->gparts[k].x[2] != s->gparts[k].part->x[2])
        error("Linked particles are not at the same position !");
    }
  }
  for (size_t k = 0; k < nr_parts; ++k) {

    if (s->parts[k].gpart != NULL) {

      if (s->parts[k].gpart->part != &s->parts[k]) error("Linking problem !");
    }
  }

586
587
  /* Hook the cells up to the parts. */
  // tic = getticks();
588
589
590
  struct part *finger = s->parts;
  struct xpart *xfinger = s->xparts;
  struct gpart *gfinger = s->gparts;
591
592
  for (int k = 0; k < s->nr_cells; k++) {
    struct cell *restrict c = &cells[k];
593
594
595
596
597
598
599
    c->parts = finger;
    c->xparts = xfinger;
    c->gparts = gfinger;
    finger = &finger[c->count];
    xfinger = &xfinger[c->count];
    gfinger = &gfinger[c->gcount];
  }
600
  // message( "hooking up cells took %.3f %s." ,
Matthieu Schaller's avatar
Matthieu Schaller committed
601
  // clocks_from_ticks(getticks() - tic), clocks_getunit());
602
603
604

  /* At this point, we have the upper-level cells, old or new. Now make
     sure that the parts in each cell are ok. */
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
  space_split(s, cells, verbose);

  if (verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
}

/**
 * @brief Split particles between cells of a hierarchy
 *
 * @param s The #space.
 * @param cells The cell hierarchy
 * @param verbose Are we talkative ?
 */
void space_split(struct space *s, struct cell *cells, int verbose) {

Matthieu Schaller's avatar
Matthieu Schaller committed
621
  const ticks tic = getticks();
622
623

  for (int k = 0; k < s->nr_cells; k++)
624
625
    scheduler_addtask(&s->e->sched, task_type_split_cell, task_subtype_none, k,
                      0, &cells[k], NULL, 0);
626
  engine_launch(s->e, s->e->nr_threads, 1 << task_type_split_cell, 0);
627

628
629
630
  if (verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
631
}
632

633
/**
634
635
 * @brief Sort the particles and condensed particles according to the given
 *indices.
636
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
637
 * @param s The #space.
638
639
640
641
 * @param ind The indices with respect to which the parts are sorted.
 * @param N The number of parts
 * @param min Lowest index.
 * @param max highest index.
642
 * @param verbose Are we talkative ?
643
 */
644

645
void space_parts_sort(struct space *s, int *ind, size_t N, int min, int max,
646
647
                      int verbose) {

Matthieu Schaller's avatar
Matthieu Schaller committed
648
  const ticks tic = getticks();
649
650

  /*Populate the global parallel_sort structure with the input data */
651
652
653
  space_sort_struct.parts = s->parts;
  space_sort_struct.xparts = s->xparts;
  space_sort_struct.ind = ind;
654
  space_sort_struct.stack_size = 2 * (max - min + 1) + 10 + s->e->nr_threads;
655
656
657
658
659
660
  if ((space_sort_struct.stack = malloc(sizeof(struct qstack) *
                                        space_sort_struct.stack_size)) == NULL)
    error("Failed to allocate sorting stack.");
  for (int i = 0; i < space_sort_struct.stack_size; i++)
    space_sort_struct.stack[i].ready = 0;

661
  /* Add the first interval. */
662
663
664
665
666
667
668
669
670
  space_sort_struct.stack[0].i = 0;
  space_sort_struct.stack[0].j = N - 1;
  space_sort_struct.stack[0].min = min;
  space_sort_struct.stack[0].max = max;
  space_sort_struct.stack[0].ready = 1;
  space_sort_struct.first = 0;
  space_sort_struct.last = 1;
  space_sort_struct.waiting = 1;

671
  /* Launch the sorting tasks. */
672
  engine_launch(s->e, s->e->nr_threads, (1 << task_type_part_sort), 0);
673
674

  /* Verify space_sort_struct. */
675
  /* for (int i = 1; i < N; i++)
676
    if (ind[i - 1] > ind[i])
677
678
      error("Sorting failed (ind[%i]=%i,ind[%i]=%i), min=%i, max=%i.", i - 1,
  ind[i - 1], i,
679
680
            ind[i], min, max);
  message("Sorting succeeded."); */
681

682
  /* Clean up. */
683
  free(space_sort_struct.stack);
684
685
686
687

  if (verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
688
}
689

690
void space_do_parts_sort() {
691

692
  /* Pointers to the sorting data. */
693
  int *ind = space_sort_struct.ind;
694
695
  struct part *parts = space_sort_struct.parts;
  struct xpart *xparts = space_sort_struct.xparts;
696

697
  /* Main loop. */
698
  while (space_sort_struct.waiting) {
699

700
    /* Grab an interval off the queue. */
701
702
    int qid =
        atomic_inc(&space_sort_struct.first) % space_sort_struct.stack_size;
703

704
    /* Wait for the entry to be ready, or for the sorting do be done. */
705
706
    while (!space_sort_struct.stack[qid].ready)
      if (!space_sort_struct.waiting) return;
707

708
    /* Get the stack entry. */
709
710
    ptrdiff_t i = space_sort_struct.stack[qid].i;
    ptrdiff_t j = space_sort_struct.stack[qid].j;
711
712
    int min = space_sort_struct.stack[qid].min;
    int max = space_sort_struct.stack[qid].max;
713
    space_sort_struct.stack[qid].ready = 0;
714

715
716
    /* Loop over sub-intervals. */
    while (1) {
717

718
      /* Bring beer. */
719
      const int pivot = (min + max) / 2;
720
721
      /* message("Working on interval [%i,%i] with min=%i, max=%i, pivot=%i.",
              i, j, min, max, pivot); */
722
723

      /* One pass of QuickSort's partitioning. */
724
725
      ptrdiff_t ii = i;
      ptrdiff_t jj = j;
726
727
728
729
      while (ii < jj) {
        while (ii <= j && ind[ii] <= pivot) ii++;
        while (jj >= i && ind[jj] > pivot) jj--;
        if (ii < jj) {
730
          size_t temp_i = ind[ii];
731
732
          ind[ii] = ind[jj];
          ind[jj] = temp_i;
733
          struct part temp_p = parts[ii];
734
735
          parts[ii] = parts[jj];
          parts[jj] = temp_p;
736
          struct xpart temp_xp = xparts[ii];
737
738
739
740
          xparts[ii] = xparts[jj];
          xparts[jj] = temp_xp;
        }
      }
741

742
      /* Verify space_sort_struct. */
743
744
745
746
747
748
749
750
751
752
753
754
      /* for (int k = i; k <= jj; k++)
        if (ind[k] > pivot) {
          message("sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%i, j=%i.", k,
                  ind[k], pivot, i, j);
          error("Partition failed (<=pivot).");
        }
      for (int k = jj + 1; k <= j; k++)
        if (ind[k] <= pivot) {
          message("sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%i, j=%i.", k,
                  ind[k], pivot, i, j);
          error("Partition failed (>pivot).");
        } */
755
756
757
758
759
760

      /* Split-off largest interval. */
      if (jj - i > j - jj + 1) {

        /* Recurse on the left? */
        if (jj > i && pivot > min) {
761
762
          qid = atomic_inc(&space_sort_struct.last) %
                space_sort_struct.stack_size;
763
764
          while (space_sort_struct.stack[qid].ready)
            ;
765
766
767
768
769
770
771
          space_sort_struct.stack[qid].i = i;
          space_sort_struct.stack[qid].j = jj;
          space_sort_struct.stack[qid].min = min;
          space_sort_struct.stack[qid].max = pivot;
          if (atomic_inc(&space_sort_struct.waiting) >=
              space_sort_struct.stack_size)
            error("Qstack overflow.");
772
          space_sort_struct.stack[qid].ready = 1;
773
        }
774

775
776
777
778
779
780
781
782
783
784
        /* Recurse on the right? */
        if (jj + 1 < j && pivot + 1 < max) {
          i = jj + 1;
          min = pivot + 1;
        } else
          break;

      } else {

        /* Recurse on the right? */
785
        if (pivot + 1 < max) {
786
787
          qid = atomic_inc(&space_sort_struct.last) %
                space_sort_struct.stack_size;
788
789
          while (space_sort_struct.stack[qid].ready)
            ;
790
791
792
793
794
795
796
          space_sort_struct.stack[qid].i = jj + 1;
          space_sort_struct.stack[qid].j = j;
          space_sort_struct.stack[qid].min = pivot + 1;
          space_sort_struct.stack[qid].max = max;
          if (atomic_inc(&space_sort_struct.waiting) >=
              space_sort_struct.stack_size)
            error("Qstack overflow.");
797
          space_sort_struct.stack[qid].ready = 1;
798
        }
799

800
801
802
803
804
805
806
        /* Recurse on the left? */
        if (jj > i && pivot > min) {
          j = jj;
          max = pivot;
        } else
          break;
      }
807

808
809
    } /* loop over sub-intervals. */

810
    atomic_dec(&space_sort_struct.waiting);
811
812

  } /* main loop. */
813
814
}

815
816
817
818
819
/**
 * @brief Sort the g-particles and condensed particles according to the given
 *indices.
 *
 * @param s The #space.
Matthieu Schaller's avatar
Matthieu Schaller committed
820
821
 * @param ind The indices with respect to which the gparts are sorted.
 * @param N The number of gparts
822
823
824
825
 * @param min Lowest index.
 * @param max highest index.
 * @param verbose Are we talkative ?
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
826
void space_gparts_sort(struct space *s, int *ind, size_t N, int min, int max,
827
828
                       int verbose) {

Matthieu Schaller's avatar
Matthieu Schaller committed
829
  const ticks tic = getticks();
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872

  /*Populate the global parallel_sort structure with the input data */
  space_sort_struct.gparts = s->gparts;
  space_sort_struct.ind = ind;
  space_sort_struct.stack_size = 2 * (max - min + 1) + 10 + s->e->nr_threads;
  if ((space_sort_struct.stack = malloc(sizeof(struct qstack) *
                                        space_sort_struct.stack_size)) == NULL)
    error("Failed to allocate sorting stack.");
  for (int i = 0; i < space_sort_struct.stack_size; i++)
    space_sort_struct.stack[i].ready = 0;

  /* Add the first interval. */
  space_sort_struct.stack[0].i = 0;
  space_sort_struct.stack[0].j = N - 1;
  space_sort_struct.stack[0].min = min;
  space_sort_struct.stack[0].max = max;
  space_sort_struct.stack[0].ready = 1;
  space_sort_struct.first = 0;
  space_sort_struct.last = 1;
  space_sort_struct.waiting = 1;

  /* Launch the sorting tasks. */
  engine_launch(s->e, s->e->nr_threads, (1 << task_type_gpart_sort), 0);

  /* Verify space_sort_struct. */
  /* for (int i = 1; i < N; i++)
    if (ind[i - 1] > ind[i])
      error("Sorting failed (ind[%i]=%i,ind[%i]=%i), min=%i, max=%i.", i - 1,
  ind[i - 1], i,
            ind[i], min, max);
  message("Sorting succeeded."); */

  /* Clean up. */
  free(space_sort_struct.stack);

  if (verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
}

void space_do_gparts_sort() {

  /* Pointers to the sorting data. */
Matthieu Schaller's avatar
Matthieu Schaller committed
873
  int *ind = space_sort_struct.ind;
874
  struct gpart *gparts = space_sort_struct.gparts;
875

876
  /* Main loop. */
877
  while (space_sort_struct.waiting) {
878

879
    /* Grab an interval off the queue. */
880
881
882
883
884
885
    int qid =
        atomic_inc(&space_sort_struct.first) % space_sort_struct.stack_size;

    /* Wait for the entry to be ready, or for the sorting do be done. */
    while (!space_sort_struct.stack[qid].ready)
      if (!space_sort_struct.waiting) return;
886

887
    /* Get the stack entry. */
888
889
890
891
892
    ptrdiff_t i = space_sort_struct.stack[qid].i;
    ptrdiff_t j = space_sort_struct.stack[qid].j;
    int min = space_sort_struct.stack[qid].min;
    int max = space_sort_struct.stack[qid].max;
    space_sort_struct.stack[qid].ready = 0;
893
894
895

    /* Loop over sub-intervals. */
    while (1) {
896

897
      /* Bring beer. */
898
899
900
      const int pivot = (min + max) / 2;
      /* message("Working on interval [%i,%i] with min=%i, max=%i, pivot=%i.",
              i, j, min, max, pivot); */
901
902

      /* One pass of QuickSort's partitioning. */
903
904
      ptrdiff_t ii = i;
      ptrdiff_t jj = j;
905
906
907
908
      while (ii < jj) {
        while (ii <= j && ind[ii] <= pivot) ii++;
        while (jj >= i && ind[jj] > pivot) jj--;
        if (ii < jj) {
909
          size_t temp_i = ind[ii];
910
911
          ind[ii] = ind[jj];
          ind[jj] = temp_i;
912
          struct gpart temp_p = gparts[ii];
913
914
915
916
          gparts[ii] = gparts[jj];
          gparts[jj] = temp_p;
        }
      }
917

918
      /* Verify space_sort_struct. */
919
920
921
922
923
924
925
926
927
928
929
930
      /* for (int k = i; k <= jj; k++)
        if (ind[k] > pivot) {
          message("sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%i, j=%i.", k,
                  ind[k], pivot, i, j);
          error("Partition failed (<=pivot).");
        }
      for (int k = jj + 1; k <= j; k++)
        if (ind[k] <= pivot) {
          message("sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%i, j=%i.", k,
                  ind[k], pivot, i, j);
          error("Partition failed (>pivot).");
        } */
931
932
933
934
935
936

      /* Split-off largest interval. */
      if (jj - i > j - jj + 1) {

        /* Recurse on the left? */
        if (jj > i && pivot > min) {
937
938
939
940
941
942
943
944
945
946
947
948
          qid = atomic_inc(&space_sort_struct.last) %
                space_sort_struct.stack_size;
          while (space_sort_struct.stack[qid].ready)
            ;
          space_sort_struct.stack[qid].i = i;
          space_sort_struct.stack[qid].j = jj;
          space_sort_struct.stack[qid].min = min;
          space_sort_struct.stack[qid].max = pivot;
          if (atomic_inc(&space_sort_struct.waiting) >=
              space_sort_struct.stack_size)
            error("Qstack overflow.");
          space_sort_struct.stack[qid].ready = 1;
949
        }
950

951
952
953
954
955
956
        /* Recurse on the right? */
        if (jj + 1 < j && pivot + 1 < max) {
          i = jj + 1;
          min = pivot + 1;
        } else
          break;
957

958
959
960
      } else {

        /* Recurse on the right? */
961
        if (pivot + 1 < max) {
962
963
964
965
966
967
968
969
970
971
972
973
          qid = atomic_inc(&space_sort_struct.last) %
                space_sort_struct.stack_size;
          while (space_sort_struct.stack[qid].ready)
            ;
          space_sort_struct.stack[qid].i = jj + 1;
          space_sort_struct.stack[qid].j = j;
          space_sort_struct.stack[qid].min = pivot + 1;
          space_sort_struct.stack[qid].max = max;
          if (atomic_inc(&space_sort_struct.waiting) >=
              space_sort_struct.stack_size)
            error("Qstack overflow.");
          space_sort_struct.stack[qid].ready = 1;
974
975
976
977
978
979
980
981
982
983
984
985
        }

        /* Recurse on the left? */
        if (jj > i && pivot > min) {
          j = jj;
          max = pivot;
        } else
          break;
      }

    } /* loop over sub-intervals. */

986
    atomic_dec(&space_sort_struct.waiting);
987
988

  } /* main loop. */
989
}
990

Pedro Gonnet's avatar
Pedro Gonnet committed
991
/**
992
 * @brief Mapping function to free the sorted indices buffers.
Pedro Gonnet's avatar
Pedro Gonnet committed
993
994
 */

995
void space_map_clearsort(struct cell *c, void *data) {
Pedro Gonnet's avatar
Pedro Gonnet committed
996

997
998
999
1000
1001
  if (c->sort != NULL) {
    free(c->sort);
    c->sort = NULL;
  }
}
Pedro Gonnet's avatar
Pedro Gonnet committed
1002

1003
1004
1005
/**
 * @brief Map a function to all particles in a cell recursively.
 *
1006
 * @param c The #cell we are working in.
1007
1008
1009
1010
 * @param fun Function pointer to apply on the cells.
 * @param data Data passed to the function fun.
 */

Pedro Gonnet's avatar
Pedro Gonnet committed
1011
1012
1013
1014
static void rec_map_parts(struct cell *c,
                          void (*fun)(struct part *p, struct cell *c,
                                      void *data),
                          void *data) {
1015
1016
1017
1018
1019
1020

  int k;

  /* No progeny? */
  if (!c->split)
    for (k = 0; k < c->count; k++) fun(&c->parts[k], c, data);
Pedro Gonnet's avatar
Pedro Gonnet committed
1021

1022
1023
1024
1025
1026
1027
  /* Otherwise, recurse. */
  else
    for (k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) rec_map_parts(c->progeny[k], fun, data);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
1028
/**
1029
 * @brief Map a function to all particles in a space.
Pedro Gonnet's avatar
Pedro Gonnet committed
1030
1031
 *
 * @param s The #space we are working in.
1032
1033
 * @param fun Function pointer to apply on the cells.
 * @param data Data passed to the function fun.
Pedro Gonnet's avatar
Pedro Gonnet committed
1034
1035
 */

1036
1037
1038
void space_map_parts(struct space *s,
                     void (*fun)(struct part *p, struct cell *c, void *data),
                     void *data) {
Pedro Gonnet's avatar
Pedro Gonnet committed
1039

1040
1041
  int cid = 0;

1042
  /* Call the recursive function on all higher-level cells. */
Pedro Gonnet's avatar
Pedro Gonnet committed
1043
1044
  for (cid = 0; cid < s->nr_cells; cid++)
    rec_map_parts(&s->cells[cid], fun, data);
1045
}
1046

1047
1048
1049
1050
1051
1052
1053
1054
/**
 * @brief Map a function to all particles in a cell recursively.
 *
 * @param c The #cell we are working in.
 * @param fun Function pointer to apply on the cells.
 */

static void rec_map_parts_xparts(struct cell *c,
1055
1056
                                 void (*fun)(struct part *p, struct xpart *xp,
                                             struct cell *c)) {
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077

  int k;

  /* No progeny? */
  if (!c->split)
    for (k = 0; k < c->count; k++) fun(&c->parts[k], &c->xparts[k], c);

  /* Otherwise, recurse. */
  else
    for (k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) rec_map_parts_xparts(c->progeny[k], fun);
}

/**
 * @brief Map a function to all particles (#part and #xpart) in a space.
 *
 * @param s The #space we are working in.
 * @param fun Function pointer to apply on the particles in the cells.
 */

void space_map_parts_xparts(struct space *s,
1078
1079
                            void (*fun)(struct part *p, struct xpart *xp,
                                        struct cell *c)) {
1080
1081
1082
1083
1084
1085
1086
1087

  int cid = 0;

  /* Call the recursive function on all higher-level cells. */
  for (cid = 0; cid < s->nr_cells; cid++)
    rec_map_parts_xparts(&s->cells[cid], fun);
}

1088
1089
1090
/**
 * @brief Map a function to all particles in a cell recursively.
 *