space.c 46.6 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 "tools.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
51

52
53
54
/* Shared sort structure. */
struct parallel_sort space_sort_struct;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

208
209
210
/* 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. */
211
#ifdef WITH_MPI
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
237
238
  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;
        }
      }
    }
  }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Peter W. Draper's avatar
Peter W. Draper committed
461
  /* Check that all parts are in the correct places. */
462
463
464
465
466
467
468
469
470
471
472
  /*  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");
    }
  }*/

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
483
484
485
486
      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
487
488
489
      const int t = gind[k];
      gind[k] = gind[nr_gparts];
      gind[nr_gparts] = t;
Matthieu Schaller's avatar
Matthieu Schaller committed
490
    } else {
491
492
493
494
      /* Increment when not exchanging otherwise we need to retest "k".*/
      k++;
    }
  }
495

496
497
498
499
500
501
502
503
504
505
506
507
508
  /* 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");
    }
  }*/

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

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

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

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

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

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

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

563
564
565
566
#ifdef WITH_MPI

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

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

587
#endif
588
589

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

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

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

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

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

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

  for (int k = 0; k < s->nr_cells; k++)
659
660
    scheduler_addtask(&s->e->sched, task_type_split_cell, task_subtype_none, k,
                      0, &cells[k], NULL, 0);
661
  engine_launch(s->e, s->e->nr_threads, 1 << task_type_split_cell, 0);
662

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

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

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

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

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

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

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

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

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

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

725
void space_do_parts_sort() {
726

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

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

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

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

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

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

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

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

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

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

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

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

      } else {

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

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

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

845
    atomic_dec(&space_sort_struct.waiting);
846
847

  } /* main loop. */
848
849
}

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

Matthieu Schaller's avatar
Matthieu Schaller committed
864
  const ticks tic = getticks();
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
907

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

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

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

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

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

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

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

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

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

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

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

993
994
995
      } else {

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

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

    } /* loop over sub-intervals. */

1021
    atomic_dec(&space_sort_struct.waiting);
1022
1023

  } /* main loop. */
1024
}
1025

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

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

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

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

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

  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
1056

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

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

1075
1076
  int cid = 0;

1077
  /* Call the recursive function on all higher-level cells. */
Pedro Gonnet's avatar
Pedro Gonnet committed
1078
1079