space.c 47.1 KB
Newer Older
1
/*******************************************************************************
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 * This file is part of SWIFT.
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
 *               2016 John A. Regan (john.a.regan@durham.ac.uk)
 *                    Tom Theuns (tom.theuns@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/>.
 *
 ******************************************************************************/
Pedro Gonnet's avatar
Pedro Gonnet committed
23
24
25
26
27
28
29
30

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

/* Some standard headers. */
#include <float.h>
#include <limits.h>
#include <math.h>
31
#include <stdlib.h>
32
#include <string.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
33

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

39
40
41
/* This object's header. */
#include "space.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
42
/* Local headers. */
43
#include "atomic.h"
44
#include "engine.h"
45
#include "error.h"
46
#include "kernel_hydro.h"
47
#include "lock.h"
48
#include "minmax.h"
49
#include "runner.h"
50
#include "threadpool.h"
51
#include "tools.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
52
53
54

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

/* Map shift vector to sortlist. */
const int sortlistID[27] = {
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
87
    /* ( -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};

88
89
90
91
92
93
94
95
96
97
98
99
/**
 * @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.
 */

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

  /* Get the relative distance between the pairs, wrapping. */
104
105
106
  const int periodic = s->periodic;
  double dx[3];
  for (int k = 0; k < 3; k++) {
107
108
109
110
111
112
113
114
115
116
117
    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. */
118
  int sid = 0;
119
  for (int k = 0; k < 3; k++)
120
121
122
123
    sid = 3 * sid + ((dx[k] < 0.0) ? 0 : ((dx[k] > 0.0) ? 2 : 1));

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

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

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

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

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

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

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

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

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

/* 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)
188
      error("Failed to aggregate the rebuild flag across nodes.");
189
190
191
    h_max = buff;
  }
#endif
192
  if (verbose) message("h_max is %.3e (cell_max=%.3e).", h_max, cell_max);
193
194

  /* Get the new putative cell dimensions. */
195
  int cdim[3];
196
  for (int k = 0; k < 3; k++)
197
198
199
200
201
202
203
204
205
    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.");

206
207
208
/* 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. */
209
#ifdef WITH_MPI
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
  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;
        }
      }
    }
  }

237
238
239
240
241
242
243
244
245
#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) {
246
      for (int k = 0; k < s->nr_cells; k++) {
247
248
249
250
251
252
253
254
        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. */
255
    for (int k = 0; k < 3; k++) {
256
257
258
259
      s->cdim[k] = cdim[k];
      s->h[k] = s->dim[k] / cdim[k];
      s->ih[k] = 1.0 / s->h[k];
    }
260
    const float dmin = fminf(s->h[0], fminf(s->h[1], s->h[2]));
261
262
263
264
265
266
267

    /* 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));
268
    for (int k = 0; k < s->nr_cells; k++)
269
270
271
      if (lock_init(&s->cells[k].lock) != 0) error("Failed to init spinlock.");

    /* Set the cell location and sizes. */
272
273
274
    for (int i = 0; i < cdim[0]; i++)
      for (int j = 0; j < cdim[1]; j++)
        for (int k = 0; k < cdim[2]; k++) {
275
276
277
278
279
280
281
282
283
284
285
286
287
          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
288
        }
289
290

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

296
#ifdef WITH_MPI
297
298
299
300
301
    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)
302
303
304
        message(
            "basic cell dimensions have increased - recalculating the "
            "global partition.");
305

306
      if (!partition_space_to_space(oldh, oldcdim, oldnodeIDs, s)) {
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324

        /* 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);
325

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

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

    /* Free the old cells, if they were allocated. */
338
    for (int k = 0; k < s->nr_cells; k++) {
339
340
341
342
343
344
345
346
347
348
349
      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
350
      s->cells[k].init = NULL;
Matthieu Schaller's avatar
Matthieu Schaller committed
351
      s->cells[k].ghost = NULL;
Matthieu Schaller's avatar
Matthieu Schaller committed
352
      s->cells[k].kick = NULL;
353
      s->cells[k].super = &s->cells[k];
354
    }
355
356
    s->maxdepth = 0;
  }
357
358
359
360

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

/**
 * @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.
368
 * @param verbose Print messages to stdout or not
369
370
 *
 */
371

372
void space_rebuild(struct space *s, double cell_max, int verbose) {
373

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

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

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

Pedro Gonnet's avatar
Pedro Gonnet committed
382
383
  size_t nr_parts = s->nr_parts;
  size_t nr_gparts = s->nr_gparts;
384
385
  struct cell *restrict cells = s->cells;

Matthieu Schaller's avatar
Matthieu Schaller committed
386
387
388
  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]};
389
390
391
392

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

410
411
412
413
414
415
416
  /* 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
417
    struct gpart *restrict gp = &s->gparts[k];
418
419
420
421
422
423
424
425
426
427
428
429
    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());

430
#ifdef WITH_MPI
431

432
  /* Move non-local parts to the end of the list. */
433
  const int local_nodeID = s->e->nodeID;
434
  for (size_t k = 0; k < nr_parts;) {
435
    if (cells[ind[k]].nodeID != local_nodeID) {
436
437
      cells[ind[k]].count -= 1;
      nr_parts -= 1;
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
438
      const struct part tp = s->parts[k];
439
440
      s->parts[k] = s->parts[nr_parts];
      s->parts[nr_parts] = tp;
441
      if (s->parts[k].gpart != NULL) {
442
        s->parts[k].gpart->id_or_neg_offset = -k;
443
444
      }
      if (s->parts[nr_parts].gpart != NULL) {
445
        s->parts[nr_parts].gpart->id_or_neg_offset = -nr_parts;
446
      }
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
447
      const struct xpart txp = s->xparts[k];
448
449
      s->xparts[k] = s->xparts[nr_parts];
      s->xparts[nr_parts] = txp;
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
450
      const int t = ind[k];
451
452
      ind[k] = ind[nr_parts];
      ind[nr_parts] = t;
Matthieu Schaller's avatar
Matthieu Schaller committed
453
    } else {
454
455
456
457
      /* Increment when not exchanging otherwise we need to retest "k".*/
      k++;
    }
  }
458

459
#ifdef SWIFT_DEBUG_CHECKS
Peter W. Draper's avatar
Peter W. Draper committed
460
  /* Check that all parts are in the correct places. */
461
  for (size_t k = 0; k < nr_parts; k++) {
462
463
464
465
466
467
468
469
    if (cells[ind[k]].nodeID != local_nodeID) {
      error("Failed to move all non-local parts to send list");
    }
  }
  for (size_t k = nr_parts; k < s->nr_parts; k++) {
    if (cells[ind[k]].nodeID == local_nodeID) {
      error("Failed to remove local parts from send list");
    }
470
471
  }
#endif
472

473
  /* Move non-local gparts to the end of the list. */
474
  for (int k = 0; k < nr_gparts;) {
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
475
476
    if (cells[gind[k]].nodeID != local_nodeID) {
      cells[gind[k]].gcount -= 1;
477
      nr_gparts -= 1;
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
478
      const struct gpart tp = s->gparts[k];
479
480
      s->gparts[k] = s->gparts[nr_gparts];
      s->gparts[nr_gparts] = tp;
481
482
      if (s->gparts[k].id_or_neg_offset <= 0) {
        s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
483
      }
484
485
486
      if (s->gparts[nr_gparts].id_or_neg_offset <= 0) {
        s->parts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
            &s->gparts[nr_gparts];
487
      }
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
488
489
490
      const int t = gind[k];
      gind[k] = gind[nr_gparts];
      gind[nr_gparts] = t;
Matthieu Schaller's avatar
Matthieu Schaller committed
491
    } else {
492
493
494
495
      /* Increment when not exchanging otherwise we need to retest "k".*/
      k++;
    }
  }
496

497
#ifdef SWIFT_DEBUG_CHECKS
498
499
500
501
502
503
504
505
506
507
  /* Check that all gparts are in the correct place (untested). */
  for (size_t k = 0; k < nr_gparts; k++) {
    if (cells[gind[k]].nodeID != local_nodeID) {
      error("Failed to move all non-local gparts to send list");
    }
  }
  for (size_t k = nr_gparts; k < s->nr_gparts; k++) {
    if (cells[gind[k]].nodeID == local_nodeID) {
      error("Failed to remove local gparts from send list");
    }
508
509
  }
#endif
510

511
512
  /* Exchange the strays, note that this potentially re-allocates
     the parts arrays. */
513
  size_t nr_parts_exchanged = s->nr_parts - nr_parts;
514
  size_t nr_gparts_exchanged = s->nr_gparts - nr_gparts;
Pedro Gonnet's avatar
Pedro Gonnet committed
515
516
517
518
  engine_exchange_strays(s->e, nr_parts, &ind[nr_parts], &nr_parts_exchanged,
                         nr_gparts, &gind[nr_gparts], &nr_gparts_exchanged);

  /* Set the new particle counts. */
519
  s->nr_parts = nr_parts + nr_parts_exchanged;
520
  s->nr_gparts = nr_gparts + nr_gparts_exchanged;
521
522

  /* Re-allocate the index array if needed.. */
523
  if (s->nr_parts > ind_size) {
524
525
    int *ind_new;
    if ((ind_new = (int *)malloc(sizeof(int) * s->nr_parts)) == NULL)
526
      error("Failed to allocate temporary particle indices.");
527
    memcpy(ind_new, ind, sizeof(int) * nr_parts);
528
529
    free(ind);
    ind = ind_new;
530
531
532
  }

  /* Assign each particle to its cell. */
Pedro Gonnet's avatar
Pedro Gonnet committed
533
  for (size_t k = nr_parts; k < s->nr_parts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
534
    const struct part *const p = &s->parts[k];
535
    ind[k] =
536
        cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
537
    cells[ind[k]].count += 1;
538
539
540
541
542
#ifdef SWIFT_DEBUG_CHECKS
    if (cells[ind[k]].nodeID != local_nodeID)
      error("Received part that does not belong to me (nodeID=%i).",
            cells[ind[k]].nodeID);
#endif
543
  }
544
  nr_parts = s->nr_parts;
545
546

#endif /* WITH_MPI */
547
548

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

  /* Re-link the gparts. */
552
  part_relink_gparts(s->parts, nr_parts, 0);
553

554
#ifdef SWIFT_DEBUG_CHECKS
555
  /* Verify space_sort_struct. */
556
557
558
559
560
561
562
563
564
565
  for (size_t k = 1; k < nr_parts; k++) {
    if (ind[k - 1] > ind[k]) {
      error("Sort failed!");
    } else if (ind[k] != cell_getid(cdim, s->parts[k].x[0] * ih[0],
                                    s->parts[k].x[1] * ih[1],
                                    s->parts[k].x[2] * ih[2])) {
      error("Incorrect indices!");
    }
  }
#endif
566
567

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

570
571
572
573
#ifdef WITH_MPI

  /* Re-allocate the index array if needed.. */
  if (s->nr_gparts > gind_size) {
574
575
    int *gind_new;
    if ((gind_new = (int *)malloc(sizeof(int) * s->nr_gparts)) == NULL)
576
      error("Failed to allocate temporary g-particle indices.");
577
    memcpy(gind_new, gind, sizeof(int) * nr_gparts);
578
579
580
581
582
    free(gind);
    gind = gind_new;
  }

  /* Assign each particle to its cell. */
583
  for (int k = nr_gparts; k < s->nr_gparts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
584
    const struct gpart *const p = &s->gparts[k];
585
586
    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
587
    cells[gind[k]].gcount += 1;
588
589
590
591
592
    /* 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;
593

594
#endif
595
596

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

  /* Re-link the parts. */
600
  part_relink_parts(s->gparts, nr_gparts, s->parts);
601
602

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

605
#ifdef SWIFT_DEBUG_CHECKS
606
607
608
  /* Verify that the links are correct */
  for (size_t k = 0; k < nr_gparts; ++k) {

609
    if (s->gparts[k].id_or_neg_offset < 0) {
610

611
      const struct part *part = &s->parts[-s->gparts[k].id_or_neg_offset];
612

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

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

622
623
624
    if (s->parts[k].gpart != NULL &&
        s->parts[k].gpart->id_or_neg_offset != -k) {
      error("Linking problem !");
625
626
    }
  }
627
#endif
628

629
630
  /* Hook the cells up to the parts. */
  // tic = getticks();
631
632
633
  struct part *finger = s->parts;
  struct xpart *xfinger = s->xparts;
  struct gpart *gfinger = s->gparts;
634
635
  for (int k = 0; k < s->nr_cells; k++) {
    struct cell *restrict c = &cells[k];
636
637
638
639
640
641
642
    c->parts = finger;
    c->xparts = xfinger;
    c->gparts = gfinger;
    finger = &finger[c->count];
    xfinger = &xfinger[c->count];
    gfinger = &gfinger[c->gcount];
  }
643
  // message( "hooking up cells took %.3f %s." ,
Matthieu Schaller's avatar
Matthieu Schaller committed
644
  // clocks_from_ticks(getticks() - tic), clocks_getunit());
645
646
647

  /* At this point, we have the upper-level cells, old or new. Now make
     sure that the parts in each cell are ok. */
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
  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
664
  const ticks tic = getticks();
665

666
  threadpool_map(&s->e->threadpool, space_split_mapper, cells, s->nr_cells,
667
                 sizeof(struct cell), 1, s);
668

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

674
/**
675
676
 * @brief Sort the particles and condensed particles according to the given
 *indices.
677
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
678
 * @param s The #space.
679
680
681
682
 * @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.
683
 * @param verbose Are we talkative ?
684
 */
685

686
void space_parts_sort(struct space *s, int *ind, size_t N, int min, int max,
687
688
                      int verbose) {

Matthieu Schaller's avatar
Matthieu Schaller committed
689
  const ticks tic = getticks();
690

691
692
693
694
695
696
  /* Populate a parallel_sort structure with the input data */
  struct parallel_sort sort_struct;
  sort_struct.parts = s->parts;
  sort_struct.xparts = s->xparts;
  sort_struct.ind = ind;
  sort_struct.stack_size = 2 * (max - min + 1) + 10 + s->e->nr_threads;
697
698
  if ((sort_struct.stack =
           malloc(sizeof(struct qstack) * sort_struct.stack_size)) == NULL)
699
    error("Failed to allocate sorting stack.");
700
701
  for (int i = 0; i < sort_struct.stack_size; i++)
    sort_struct.stack[i].ready = 0;
702

703
  /* Add the first interval. */
704
705
706
707
708
709
710
711
712
713
714
  sort_struct.stack[0].i = 0;
  sort_struct.stack[0].j = N - 1;
  sort_struct.stack[0].min = min;
  sort_struct.stack[0].max = max;
  sort_struct.stack[0].ready = 1;
  sort_struct.first = 0;
  sort_struct.last = 1;
  sort_struct.waiting = 1;

  /* Launch the sorting tasks with a stride of zero such that the same
     map data is passed to each thread. */
715
716
  threadpool_map(&s->e->threadpool, space_parts_sort_mapper, &sort_struct,
                 s->e->threadpool.num_threads, 0, 1, NULL);
717

718
#ifdef SWIFT_DEBUG_CHECKS
719
  /* Verify space_sort_struct. */
720
  for (int i = 1; i < N; i++)
721
    if (ind[i - 1] > ind[i])
722
      error("Sorting failed (ind[%i]=%i,ind[%i]=%i), min=%i, max=%i.", i - 1,
723
724
725
            ind[i - 1], i, ind[i], min, max);
  message("Sorting succeeded.");
#endif
726

727
  /* Clean up. */
728
  free(sort_struct.stack);
729
730
731
732

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

735
736
void space_parts_sort_mapper(void *map_data, int num_elements,
                             void *extra_data) {
737
738
739

  /* Unpack the mapping data. */
  struct parallel_sort *sort_struct = (struct parallel_sort *)map_data;
740

741
  /* Pointers to the sorting data. */
742
743
744
  int *ind = sort_struct->ind;
  struct part *parts = sort_struct->parts;
  struct xpart *xparts = sort_struct->xparts;
745

746
  /* Main loop. */
747
  while (sort_struct->waiting) {
748

749
    /* Grab an interval off the queue. */
750
    int qid = atomic_inc(&sort_struct->first) % sort_struct->stack_size;
751

752
    /* Wait for the entry to be ready, or for the sorting do be done. */
753
754
    while (!sort_struct->stack[qid].ready)
      if (!sort_struct->waiting) return;
755

756
    /* Get the stack entry. */
757
758
759
760
761
    ptrdiff_t i = sort_struct->stack[qid].i;
    ptrdiff_t j = sort_struct->stack[qid].j;
    int min = sort_struct->stack[qid].min;
    int max = sort_struct->stack[qid].max;
    sort_struct->stack[qid].ready = 0;
762

763
764
    /* Loop over sub-intervals. */
    while (1) {
765

766
      /* Bring beer. */
767
      const int pivot = (min + max) / 2;
768
769
      /* message("Working on interval [%i,%i] with min=%i, max=%i, pivot=%i.",
              i, j, min, max, pivot); */
770
771

      /* One pass of QuickSort's partitioning. */
772
773
      ptrdiff_t ii = i;
      ptrdiff_t jj = j;
774
775
776
777
      while (ii < jj) {
        while (ii <= j && ind[ii] <= pivot) ii++;
        while (jj >= i && ind[jj] > pivot) jj--;
        if (ii < jj) {
778
          size_t temp_i = ind[ii];
779
780
          ind[ii] = ind[jj];
          ind[jj] = temp_i;
781
          struct part temp_p = parts[ii];
782
783
          parts[ii] = parts[jj];
          parts[jj] = temp_p;
784
          struct xpart temp_xp = xparts[ii];
785
786
787
788
          xparts[ii] = xparts[jj];
          xparts[jj] = temp_xp;
        }
      }
789

790
#ifdef SWIFT_DEBUG_CHECKS
791
      /* Verify space_sort_struct. */
792
      for (int k = i; k <= jj; k++)
793
        if (ind[k] > pivot) {
794
795
          message("sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%li, j=%li.",
                  k, ind[k], pivot, i, j);
796
797
798
799
          error("Partition failed (<=pivot).");
        }
      for (int k = jj + 1; k <= j; k++)
        if (ind[k] <= pivot) {
800
801
          message("sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%li, j=%li.",
                  k, ind[k], pivot, i, j);
802
          error("Partition failed (>pivot).");
803
804
        }
#endif
805
806
807
808
809
810

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

        /* Recurse on the left? */
        if (jj > i && pivot > min) {
811
          qid = atomic_inc(&sort_struct->last) % sort_struct->stack_size;
812
          while (sort_struct->stack[qid].ready)
813
            ;
814
815
816
817
          sort_struct->stack[qid].i = i;
          sort_struct->stack[qid].j = jj;
          sort_struct->stack[qid].min = min;
          sort_struct->stack[qid].max = pivot;
818
          if (atomic_inc(&sort_struct->waiting) >= sort_struct->stack_size)
819
            error("Qstack overflow.");
820
          sort_struct->stack[qid].ready = 1;
821
        }
822

823
824
825
826
827
828
829
830
831
832
        /* Recurse on the right? */
        if (jj + 1 < j && pivot + 1 < max) {
          i = jj + 1;
          min = pivot + 1;
        } else
          break;

      } else {

        /* Recurse on the right? */
833
        if (pivot + 1 < max) {
834
          qid = atomic_inc(&sort_struct->last) % sort_struct->stack_size;
835
          while (sort_struct->stack[qid].ready)
836
            ;
837
838
839
840
          sort_struct->stack[qid].i = jj + 1;
          sort_struct->stack[qid].j = j;
          sort_struct->stack[qid].min = pivot + 1;
          sort_struct->stack[qid].max = max;
841
          if (atomic_inc(&sort_struct->waiting) >= sort_struct->stack_size)
842
            error("Qstack overflow.");
843
          sort_struct->stack[qid].ready = 1;
844
        }
845

846
847
848
849
850
851
852
        /* Recurse on the left? */
        if (jj > i && pivot > min) {
          j = jj;
          max = pivot;
        } else
          break;
      }
853

854
855
    } /* loop over sub-intervals. */

856
    atomic_dec(&sort_struct->waiting);
857
858

  } /* main loop. */
859
860
}

861
862
863
864
865
/**
 * @brief Sort the g-particles and condensed particles according to the given
 *indices.
 *
 * @param s The #space.
Matthieu Schaller's avatar
Matthieu Schaller committed
866
867
 * @param ind The indices with respect to which the gparts are sorted.
 * @param N The number of gparts
868
869
870
871
 * @param min Lowest index.
 * @param max highest index.
 * @param verbose Are we talkative ?
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
872
void space_gparts_sort(struct space *s, int *ind, size_t N, int min, int max,
873
874
                       int verbose) {

Matthieu Schaller's avatar
Matthieu Schaller committed
875
  const ticks tic = getticks();
876

877
878
879
880
881
  /*Populate a global parallel_sort structure with the input data */
  struct parallel_sort sort_struct;
  sort_struct.gparts = s->gparts;
  sort_struct.ind = ind;
  sort_struct.stack_size = 2 * (max - min + 1) + 10 + s->e->nr_threads;
882
883
  if ((sort_struct.stack =
           malloc(sizeof(struct qstack) * sort_struct.stack_size)) == NULL)
884
    error("Failed to allocate sorting stack.");
885
886
  for (int i = 0; i < sort_struct.stack_size; i++)
    sort_struct.stack[i].ready = 0;
887
888

  /* Add the first interval. */
889
890
891
892
893
894
895
896
897
898
899
  sort_struct.stack[0].i = 0;
  sort_struct.stack[0].j = N - 1;
  sort_struct.stack[0].min = min;
  sort_struct.stack[0].max = max;
  sort_struct.stack[0].ready = 1;
  sort_struct.first = 0;
  sort_struct.last = 1;
  sort_struct.waiting = 1;

  /* Launch the sorting tasks with a stride of zero such that the same
     map data is passed to each thread. */
900
901
  threadpool_map(&s->e->threadpool, space_gparts_sort_mapper, &sort_struct,
                 s->e->threadpool.num_threads, 0, 1, NULL);
902

903
#ifdef SWIFT_DEBUG_CHECKS
904
  /* Verify space_sort_struct. */
905
  for (int i = 1; i < N; i++)
906
907
    if (ind[i - 1] > ind[i])
      error("Sorting failed (ind[%i]=%i,ind[%i]=%i), min=%i, max=%i.", i - 1,
908
909
910
            ind[i - 1], i, ind[i], min, max);
  message("Sorting succeeded.");
#endif
911
912

  /* Clean up. */
913
  free(sort_struct.stack);
914
915
916
917
918
919

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

920
921
void space_gparts_sort_mapper(void *map_data, int num_elements,
                              void *extra_data) {
922
923
924

  /* Unpack the mapping data. */
  struct parallel_sort *sort_struct = (struct parallel_sort *)map_data;
925
926

  /* Pointers to the sorting data. */
927
928
  int *ind = sort_struct->ind;
  struct gpart *gparts = sort_struct->gparts;
929

930
  /* Main loop. */
931
  while (sort_struct->waiting) {
932

933
    /* Grab an interval off the queue. */
934
    int qid = atomic_inc(&sort_struct->first) % sort_struct->stack_size;
935
936

    /* Wait for the entry to be ready, or for the sorting do be done. */
937
938
    while (!sort_struct->stack[qid].ready)
      if (!sort_struct->waiting) return;
939

940
    /* Get the stack entry. */
941
942
943
944
945
    ptrdiff_t i = sort_struct->stack[qid].i;
    ptrdiff_t j = sort_struct->stack[qid].j;
    int min = sort_struct->stack[qid].min;
    int max = sort_struct->stack[qid].max;
    sort_struct->stack[qid].ready = 0;
946
947
948

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

950
      /* Bring beer. */
951
952
953
      const int pivot = (min + max) / 2;
      /* message("Working on interval [%i,%i] with min=%i, max=%i, pivot=%i.",
              i, j, min, max, pivot); */
954
955

      /* One pass of QuickSort's partitioning. */
956
957
      ptrdiff_t ii = i;
      ptrdiff_t jj = j;
958
959
960
961
      while (ii < jj) {
        while (ii <= j && ind[ii] <= pivot) ii++;
        while (jj >= i && ind[jj] > pivot) jj--;
        if (ii < jj) {
962
          size_t temp_i = ind[ii];
963
964
          ind[ii] = ind[jj];
          ind[jj] = temp_i;
965
          struct gpart temp_p = gparts[ii];
966
967
968
969
          gparts[ii] = gparts[jj];
          gparts[jj] = temp_p;
        }
      }