space.c 46.4 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 <string.h>
32
#include <stdlib.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
55
/* Shared sort structure. */
struct parallel_sort space_sort_struct;

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

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

91
92
93
94
95
96
97
98
99
100
101
102
/**
 * @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.
 */

103
104
105
106
int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
                 double *shift) {

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

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

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

138
/**
139
 * @brief Recursively dismantle a cell tree.
140
141
 *
 */
142
143
144
145

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

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

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

162
void space_regrid(struct space *s, double cell_max, int verbose) {
163

164
  const size_t nr_parts = s->nr_parts;
165
  struct cell *restrict c;
166
  ticks tic = getticks();
167
168
169

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  if (verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
365
}
366
367
368
369
370
371

/**
 * @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.
372
 * @param verbose Print messages to stdout or not
373
374
 *
 */
375

376
void space_rebuild(struct space *s, double cell_max, int verbose) {
377

Matthieu Schaller's avatar
Matthieu Schaller committed
378
  const ticks tic = getticks();
379
380

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

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

Pedro Gonnet's avatar
Pedro Gonnet committed
386
387
  size_t nr_parts = s->nr_parts;
  size_t nr_gparts = s->nr_gparts;
388
389
  struct cell *restrict cells = s->cells;

Matthieu Schaller's avatar
Matthieu Schaller committed
390
391
392
  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]};
393
394
395
396

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

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

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

Peter W. Draper's avatar
Peter W. Draper committed
462
  /* Check that all parts are in the correct places. */
463
464
465
466
467
468
469
470
471
472
473
  /*  for (size_t k = 0; k < nr_parts; k++) {
    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");
    }
  }*/

474
  /* Move non-local gparts to the end of the list. */
475
  for (int k = 0; k < nr_gparts;) {
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
476
477
    if (cells[gind[k]].nodeID != local_nodeID) {
      cells[gind[k]].gcount -= 1;
478
      nr_gparts -= 1;
Matthieu Schaller's avatar
Bug fix    
Matthieu Schaller committed
479
      const struct gpart tp = s->gparts[k];
480
481
      s->gparts[k] = s->gparts[nr_gparts];
      s->gparts[nr_gparts] = tp;
482
483
484
485
486
487
      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
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
498
499
500
501
502
503
504
505
506
507
508
509
  /* 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");
    }
  }*/

510
511
  /* Exchange the strays, note that this potentially re-allocates
     the parts arrays. */
512
  size_t nr_parts_exchanged = s->nr_parts - nr_parts;
513
  size_t nr_gparts_exchanged = s->nr_gparts - nr_gparts;
Pedro Gonnet's avatar
Pedro Gonnet committed
514
515
516
517
  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. */
518
  s->nr_parts = nr_parts + nr_parts_exchanged;
519
  s->nr_gparts = nr_gparts + nr_gparts_exchanged;
520
521

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

  /* Assign each particle to its cell. */
Pedro Gonnet's avatar
Pedro Gonnet committed
532
  for (size_t k = nr_parts; k < s->nr_parts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
533
    const struct part *const p = &s->parts[k];
534
    ind[k] =
535
        cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
536
537
538
539
    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 ); */
540
  }
541
  nr_parts = s->nr_parts;
542
543
544
#endif

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

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

551
  /* Verify space_sort_struct. */
552
  /* for ( k = 1 ; k < nr_parts ; k++ ) {
553
      if ( ind[k-1] > ind[k] ) {
554
555
          error( "Sort failed!" );
          }
556
      else if ( ind[k] != cell_getid( cdim , parts[k].x[0]*ih[0] ,
557
558
559
560
561
     parts[k].x[1]*ih[1] , parts[k].x[2]*ih[2] ) )
          error( "Incorrect indices!" );
      } */

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

564
565
566
567
#ifdef WITH_MPI

  /* Re-allocate the index array if needed.. */
  if (s->nr_gparts > gind_size) {
568
569
    int *gind_new;
    if ((gind_new = (int *)malloc(sizeof(int) * s->nr_gparts)) == NULL)
570
      error("Failed to allocate temporary g-particle indices.");
571
    memcpy(gind_new, gind, sizeof(int) * nr_gparts);
572
573
574
575
576
    free(gind);
    gind = gind_new;
  }

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

588
#endif
589
590

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

  /* Re-link the parts. */
594
  for (int k = 0; k < nr_gparts; k++)
595
    if (s->gparts[k].id > 0) s->gparts[k].part->gpart = &s->gparts[k];
596
597

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

600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
  /* 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 !");
    }
  }

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

  /* At this point, we have the upper-level cells, old or new. Now make
     sure that the parts in each cell are ok. */
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
  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
657
  const ticks tic = getticks();
658

659
660
  threadpool_map(&s->e->threadpool, space_split_mapper, cells, s->nr_cells,
                 sizeof(struct cell), s);
661

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

667
/**
668
669
 * @brief Sort the particles and condensed particles according to the given
 *indices.
670
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
671
 * @param s The #space.
672
673
674
675
 * @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.
676
 * @param verbose Are we talkative ?
677
 */
678

679
void space_parts_sort(struct space *s, int *ind, size_t N, int min, int max,
680
681
                      int verbose) {

Matthieu Schaller's avatar
Matthieu Schaller committed
682
  const ticks tic = getticks();
683
684

  /*Populate the global parallel_sort structure with the input data */
685
686
687
  space_sort_struct.parts = s->parts;
  space_sort_struct.xparts = s->xparts;
  space_sort_struct.ind = ind;
688
  space_sort_struct.stack_size = 2 * (max - min + 1) + 10 + s->e->nr_threads;
689
690
691
692
693
694
  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;

695
  /* Add the first interval. */
696
697
698
699
700
701
702
703
704
  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;

705
  /* Launch the sorting tasks. */
706
  engine_launch(s->e, s->e->nr_threads, (1 << task_type_part_sort), 0);
707
708

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

716
  /* Clean up. */
717
  free(space_sort_struct.stack);
718
719
720
721

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

724
void space_do_parts_sort() {
725

726
  /* Pointers to the sorting data. */
727
  int *ind = space_sort_struct.ind;
728
729
  struct part *parts = space_sort_struct.parts;
  struct xpart *xparts = space_sort_struct.xparts;
730

731
  /* Main loop. */
732
  while (space_sort_struct.waiting) {
733

734
    /* Grab an interval off the queue. */
735
736
    int qid =
        atomic_inc(&space_sort_struct.first) % space_sort_struct.stack_size;
737

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

742
    /* Get the stack entry. */
743
744
    ptrdiff_t i = space_sort_struct.stack[qid].i;
    ptrdiff_t j = space_sort_struct.stack[qid].j;
745
746
    int min = space_sort_struct.stack[qid].min;
    int max = space_sort_struct.stack[qid].max;
747
    space_sort_struct.stack[qid].ready = 0;
748

749
750
    /* Loop over sub-intervals. */
    while (1) {
751

752
      /* Bring beer. */
753
      const int pivot = (min + max) / 2;
754
755
      /* message("Working on interval [%i,%i] with min=%i, max=%i, pivot=%i.",
              i, j, min, max, pivot); */
756
757

      /* One pass of QuickSort's partitioning. */
758
759
      ptrdiff_t ii = i;
      ptrdiff_t jj = j;
760
761
762
763
      while (ii < jj) {
        while (ii <= j && ind[ii] <= pivot) ii++;
        while (jj >= i && ind[jj] > pivot) jj--;
        if (ii < jj) {
764
          size_t temp_i = ind[ii];
765
766
          ind[ii] = ind[jj];
          ind[jj] = temp_i;
767
          struct part temp_p = parts[ii];
768
769
          parts[ii] = parts[jj];
          parts[jj] = temp_p;
770
          struct xpart temp_xp = xparts[ii];
771
772
773
774
          xparts[ii] = xparts[jj];
          xparts[jj] = temp_xp;
        }
      }
775

776
      /* Verify space_sort_struct. */
777
778
779
780
781
782
783
784
785
786
787
788
      /* 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).");
        } */
789
790
791
792
793
794

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

        /* Recurse on the left? */
        if (jj > i && pivot > min) {
795
796
          qid = atomic_inc(&space_sort_struct.last) %
                space_sort_struct.stack_size;
797
798
          while (space_sort_struct.stack[qid].ready)
            ;
799
800
801
802
803
804
805
          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.");
806
          space_sort_struct.stack[qid].ready = 1;
807
        }
808

809
810
811
812
813
814
815
816
817
818
        /* Recurse on the right? */
        if (jj + 1 < j && pivot + 1 < max) {
          i = jj + 1;
          min = pivot + 1;
        } else
          break;

      } else {

        /* Recurse on the right? */
819
        if (pivot + 1 < max) {
820
821
          qid = atomic_inc(&space_sort_struct.last) %
                space_sort_struct.stack_size;
822
823
          while (space_sort_struct.stack[qid].ready)
            ;
824
825
826
827
828
829
830
          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.");
831
          space_sort_struct.stack[qid].ready = 1;
832
        }
833

834
835
836
837
838
839
840
        /* Recurse on the left? */
        if (jj > i && pivot > min) {
          j = jj;
          max = pivot;
        } else
          break;
      }
841

842
843
    } /* loop over sub-intervals. */

844
    atomic_dec(&space_sort_struct.waiting);
845
846

  } /* main loop. */
847
848
}

849
850
851
852
853
/**
 * @brief Sort the g-particles and condensed particles according to the given
 *indices.
 *
 * @param s The #space.
Matthieu Schaller's avatar
Matthieu Schaller committed
854
855
 * @param ind The indices with respect to which the gparts are sorted.
 * @param N The number of gparts
856
857
858
859
 * @param min Lowest index.
 * @param max highest index.
 * @param verbose Are we talkative ?
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
860
void space_gparts_sort(struct space *s, int *ind, size_t N, int min, int max,
861
862
                       int verbose) {

Matthieu Schaller's avatar
Matthieu Schaller committed
863
  const ticks tic = getticks();
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906

  /*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
907
  int *ind = space_sort_struct.ind;
908
  struct gpart *gparts = space_sort_struct.gparts;
909

910
  /* Main loop. */
911
  while (space_sort_struct.waiting) {
912

913
    /* Grab an interval off the queue. */
914
915
916
917
918
919
    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;
920

921
    /* Get the stack entry. */
922
923
924
925
926
    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;
927
928
929

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

931
      /* Bring beer. */
932
933
934
      const int pivot = (min + max) / 2;
      /* message("Working on interval [%i,%i] with min=%i, max=%i, pivot=%i.",
              i, j, min, max, pivot); */
935
936

      /* One pass of QuickSort's partitioning. */
937
938
      ptrdiff_t ii = i;
      ptrdiff_t jj = j;
939
940
941
942
      while (ii < jj) {
        while (ii <= j && ind[ii] <= pivot) ii++;
        while (jj >= i && ind[jj] > pivot) jj--;
        if (ii < jj) {
943
          size_t temp_i = ind[ii];
944
945
          ind[ii] = ind[jj];
          ind[jj] = temp_i;
946
          struct gpart temp_p = gparts[ii];
947
948
949
950
          gparts[ii] = gparts[jj];
          gparts[jj] = temp_p;
        }
      }
951

952
      /* Verify space_sort_struct. */
953
954
955
956
957
958
959
960
961
962
963
964
      /* 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).");
        } */
965
966
967
968
969
970

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

        /* Recurse on the left? */
        if (jj > i && pivot > min) {
971
972
973
974
975
976
977
978
979
980
981
982
          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;
983
        }
984

985
986
987
988
989
990
        /* Recurse on the right? */
        if (jj + 1 < j && pivot + 1 < max) {
          i = jj + 1;
          min = pivot + 1;
        } else
          break;
991

992
993
994
      } else {

        /* Recurse on the right? */
995
        if (pivot + 1 < max) {
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
          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;
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
        }

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

    } /* loop over sub-intervals. */

1020
    atomic_dec(&space_sort_struct.waiting);
1021
1022

  } /* main loop. */
1023
}
1024

Pedro Gonnet's avatar
Pedro Gonnet committed
1025
/**
1026
 * @brief Mapping function to free the sorted indices buffers.
Pedro Gonnet's avatar
Pedro Gonnet committed
1027
1028
 */

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

1031
1032
1033
1034
1035
  if (c->sort != NULL) {
    free(c->sort);
    c->sort = NULL;
  }
}
Pedro Gonnet's avatar
Pedro Gonnet committed
1036

1037
1038
1039
/**
 * @brief Map a function to all particles in a cell recursively.
 *
1040
 * @param c The #cell we are working in.
1041
1042
1043
1044
 * @param fun Function pointer to apply on the cells.
 * @param data Data passed to the function fun.
 */

Pedro Gonnet's avatar
Pedro Gonnet committed
1045
1046
1047
1048
static void rec_map_parts(struct cell *c,
                          void (*fun)(struct part *p, struct cell *c,
                                      void *data),
                          void *data) {
1049
1050
1051
1052
1053
1054

  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
1055

1056
1057
1058
1059
1060
1061
  /* 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
1062
/**
1063
 * @brief Map a function to all particles in a space.
Pedro Gonnet's avatar
Pedro Gonnet committed
1064
1065
 *
 * @param s The #space we are working in.
1066
1067
 * @param fun Function pointer to apply on the cells.
 * @param data Data passed to the function fun.
Pedro Gonnet's avatar
Pedro Gonnet committed
1068
1069
 */

1070
1071
1072
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
1073

1074
1075
  int cid = 0;

1076
  /* Call the recursive function on all higher-level cells. */
Pedro Gonnet's avatar
Pedro Gonnet committed
1077
1078
  for (cid = 0; cid < s->nr_cells