space.c 46.9 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
211
212
213
  /* Check if we have enough cells for gravity. */
  if (s->gravity && (cdim[0] < 8 || cdim[1] < 8 || cdim[2] < 8))
    error(
        "Must have at least 8 cells in each spatial dimension when gravity "
        "is switched on.");

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

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

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

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

    /* Be verbose about the change. */
299
300
301
    if (verbose)
      message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1],
              cdim[2]);
302
303
    fflush(stdout);

304
#ifdef WITH_MPI
305
306
307
308
309
    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)
310
311
312
        message(
            "basic cell dimensions have increased - recalculating the "
            "global partition.");
313

314
      if (!partition_space_to_space(oldh, oldcdim, oldnodeIDs, s)) {
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332

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

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

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

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

  if (verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
370
}
371
372
373
374
375
376

/**
 * @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.
377
 * @param verbose Print messages to stdout or not
378
379
 *
 */
380

381
void space_rebuild(struct space *s, double cell_max, int verbose) {
382

Matthieu Schaller's avatar
Matthieu Schaller committed
383
  const ticks tic = getticks();
384
385

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

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

Pedro Gonnet's avatar
Pedro Gonnet committed
391
392
  size_t nr_parts = s->nr_parts;
  size_t nr_gparts = s->nr_gparts;
393
394
  struct cell *restrict cells = s->cells;

Matthieu Schaller's avatar
Matthieu Schaller committed
395
396
397
  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]};
398
399
400
401

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

419
420
421
422
423
424
425
  /* 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
426
    struct gpart *restrict gp = &s->gparts[k];
427
428
429
430
431
432
433
434
435
436
437
438
    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());

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

Peter W. Draper's avatar
Peter W. Draper committed
467
  /* Check that all parts are in the correct places. */
468
469
470
471
472
473
474
475
476
477
478
  /*  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");
    }
  }*/

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

502
503
504
505
506
507
508
509
510
511
512
513
514
  /* 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");
    }
  }*/

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

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

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

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

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

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

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

569
570
571
572
#ifdef WITH_MPI

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

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

593
#endif
594
595

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

  /* Re-link the parts. */
599
  for (int k = 0; k < nr_gparts; k++)
600
    if (s->gparts[k].id > 0) s->gparts[k].part->gpart = &s->gparts[k];
601
602

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

605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
  /* 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 !");
    }
  }

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

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

  for (int k = 0; k < s->nr_cells; k++)
665
666
    scheduler_addtask(&s->e->sched, task_type_split_cell, task_subtype_none, k,
                      0, &cells[k], NULL, 0);
667
  engine_launch(s->e, s->e->nr_threads, 1 << task_type_split_cell, 0);
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

  /*Populate the global parallel_sort structure with the input data */
692
693
694
  space_sort_struct.parts = s->parts;
  space_sort_struct.xparts = s->xparts;
  space_sort_struct.ind = ind;
695
  space_sort_struct.stack_size = 2 * (max - min + 1) + 10 + s->e->nr_threads;
696
697
698
699
700
701
  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;

702
  /* Add the first interval. */
703
704
705
706
707
708
709
710
711
  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;

712
  /* Launch the sorting tasks. */
713
  engine_launch(s->e, s->e->nr_threads, (1 << task_type_part_sort), 0);
714
715

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

723
  /* Clean up. */
724
  free(space_sort_struct.stack);
725
726
727
728

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

731
void space_do_parts_sort() {
732

733
  /* Pointers to the sorting data. */
734
  int *ind = space_sort_struct.ind;
735
736
  struct part *parts = space_sort_struct.parts;
  struct xpart *xparts = space_sort_struct.xparts;
737

738
  /* Main loop. */
739
  while (space_sort_struct.waiting) {
740

741
    /* Grab an interval off the queue. */
742
743
    int qid =
        atomic_inc(&space_sort_struct.first) % space_sort_struct.stack_size;
744

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

749
    /* Get the stack entry. */
750
751
    ptrdiff_t i = space_sort_struct.stack[qid].i;
    ptrdiff_t j = space_sort_struct.stack[qid].j;
752
753
    int min = space_sort_struct.stack[qid].min;
    int max = space_sort_struct.stack[qid].max;
754
    space_sort_struct.stack[qid].ready = 0;
755

756
757
    /* Loop over sub-intervals. */
    while (1) {
758

759
      /* Bring beer. */
760
      const int pivot = (min + max) / 2;
761
762
      /* message("Working on interval [%i,%i] with min=%i, max=%i, pivot=%i.",
              i, j, min, max, pivot); */
763
764

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

783
      /* Verify space_sort_struct. */
784
785
786
787
788
789
790
791
792
793
794
795
      /* 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).");
        } */
796
797
798
799
800
801

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

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

816
817
818
819
820
821
822
823
824
825
        /* Recurse on the right? */
        if (jj + 1 < j && pivot + 1 < max) {
          i = jj + 1;
          min = pivot + 1;
        } else
          break;

      } else {

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

841
842
843
844
845
846
847
        /* Recurse on the left? */
        if (jj > i && pivot > min) {
          j = jj;
          max = pivot;
        } else
          break;
      }
848

849
850
    } /* loop over sub-intervals. */

851
    atomic_dec(&space_sort_struct.waiting);
852
853

  } /* main loop. */
854
855
}

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

Matthieu Schaller's avatar
Matthieu Schaller committed
870
  const ticks tic = getticks();
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
908
909
910
911
912
913

  /*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
914
  int *ind = space_sort_struct.ind;
915
  struct gpart *gparts = space_sort_struct.gparts;
916

917
  /* Main loop. */
918
  while (space_sort_struct.waiting) {
919

920
    /* Grab an interval off the queue. */
921
922
923
924
925
926
    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;
927

928
    /* Get the stack entry. */
929
930
931
932
933
    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;
934
935
936

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

938
      /* Bring beer. */
939
940
941
      const int pivot = (min + max) / 2;
      /* message("Working on interval [%i,%i] with min=%i, max=%i, pivot=%i.",
              i, j, min, max, pivot); */
942
943

      /* One pass of QuickSort's partitioning. */
944
945
      ptrdiff_t ii = i;
      ptrdiff_t jj = j;
946
947
948
949
      while (ii < jj) {
        while (ii <= j && ind[ii] <= pivot) ii++;
        while (jj >= i && ind[jj] > pivot) jj--;
        if (ii < jj) {
950
          size_t temp_i = ind[ii];
951
952
          ind[ii] = ind[jj];
          ind[jj] = temp_i;
953
          struct gpart temp_p = gparts[ii];
954
955
956
957
          gparts[ii] = gparts[jj];
          gparts[jj] = temp_p;
        }
      }
958

959
      /* Verify space_sort_struct. */
960
961
962
963
964
965
966
967
968
969
970
971
      /* 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).");
        } */
972
973
974
975
976
977

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

        /* Recurse on the left? */
        if (jj > i && pivot > min) {
978
979
980
981
982
983
984
985
986
987
988
989
          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;
990
        }
991

992
993
994
995
996
997
        /* Recurse on the right? */
        if (jj + 1 < j && pivot + 1 < max) {
          i = jj + 1;
          min = pivot + 1;
        } else
          break;
998

999
1000
      } else {

For faster browsing, not all history is shown. View entire blame