cell.c 33.6 KB
Newer Older
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
5
6
7
 *                    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)
8
 *
9
10
11
12
 * 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.
13
 *
14
15
16
17
 * 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.
18
 *
19
20
 * 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/>.
21
 *
22
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
32
33
34
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
35

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

41
42
/* Switch off timers. */
#ifdef TIMER
43
#undef TIMER
44
45
#endif

46
47
48
/* This object's header. */
#include "cell.h"

49
/* Local headers. */
50
#include "active.h"
51
#include "atomic.h"
52
#include "drift.h"
53
#include "error.h"
54
#include "gravity.h"
55
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
56
#include "hydro_properties.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
57
#include "memswap.h"
58
#include "minmax.h"
59
#include "scheduler.h"
60
61
#include "space.h"
#include "timers.h"
62

63
64
int counter = 0;

65
66
67
/* Global variables. */
int cell_next_tag = 0;

68
69
70
71
72
/**
 * @brief Get the size of the cell subtree.
 *
 * @param c The #cell.
 */
73
int cell_getsize(struct cell *c) {
74

Pedro Gonnet's avatar
Pedro Gonnet committed
75
76
  /* Number of cells in this subtree. */
  int count = 1;
77

78
79
  /* Sum up the progeny if split. */
  if (c->split)
Pedro Gonnet's avatar
Pedro Gonnet committed
80
    for (int k = 0; k < 8; k++)
81
82
83
84
85
86
87
      if (c->progeny[k] != NULL) count += cell_getsize(c->progeny[k]);

  /* Return the final count. */
  return count;
}

/**
88
89
90
91
92
93
94
95
 * @brief Unpack the data of a given cell and its sub-cells.
 *
 * @param pc An array of packed #pcell.
 * @param c The #cell in which to unpack the #pcell.
 * @param s The #space in which the cells are created.
 *
 * @return The number of cells created.
 */
96
97
int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {

98
99
#ifdef WITH_MPI

100
101
  /* Unpack the current pcell. */
  c->h_max = pc->h_max;
102
103
  c->ti_end_min = pc->ti_end_min;
  c->ti_end_max = pc->ti_end_max;
104
  c->ti_old = pc->ti_old;
105
  c->count = pc->count;
106
  c->gcount = pc->gcount;
107
  c->tag = pc->tag;
Matthieu Schaller's avatar
Matthieu Schaller committed
108

109
110
  /* Number of new cells created. */
  int count = 1;
111
112

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
113
  for (int k = 0; k < 8; k++)
114
    if (pc->progeny[k] >= 0) {
115
116
      struct cell *temp;
      space_getcells(s, 1, &temp);
117
      temp->count = 0;
118
      temp->gcount = 0;
119
120
121
      temp->loc[0] = c->loc[0];
      temp->loc[1] = c->loc[1];
      temp->loc[2] = c->loc[2];
122
123
124
      temp->width[0] = c->width[0] / 2;
      temp->width[1] = c->width[1] / 2;
      temp->width[2] = c->width[2] / 2;
125
      temp->dmin = c->dmin / 2;
126
127
128
      if (k & 4) temp->loc[0] += temp->width[0];
      if (k & 2) temp->loc[1] += temp->width[1];
      if (k & 1) temp->loc[2] += temp->width[2];
129
130
      temp->depth = c->depth + 1;
      temp->split = 0;
131
      temp->dx_max = 0.f;
132
133
134
135
136
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
      count += cell_unpack(&pc[pc->progeny[k]], temp, s);
137
138
    }

139
  /* Return the total number of unpacked cells. */
140
  c->pcell_size = count;
141
  return count;
142
143
144
145
146

#else
  error("SWIFT was not compiled with MPI support.");
  return 0;
#endif
147
}
148

149
/**
150
 * @brief Link the cells recursively to the given #part array.
151
152
153
154
155
156
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */
157
int cell_link_parts(struct cell *c, struct part *parts) {
158

159
160
161
  c->parts = parts;

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
162
163
164
165
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
166
        offset += cell_link_parts(c->progeny[k], &parts[offset]);
Pedro Gonnet's avatar
Pedro Gonnet committed
167
168
    }
  }
169

170
  /* Return the total number of linked particles. */
171
172
  return c->count;
}
173

174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
/**
 * @brief Link the cells recursively to the given #gpart array.
 *
 * @param c The #cell.
 * @param gparts The #gpart array.
 *
 * @return The number of particles linked.
 */
int cell_link_gparts(struct cell *c, struct gpart *gparts) {

  c->gparts = gparts;

  /* Fill the progeny recursively, depth-first. */
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
        offset += cell_link_gparts(c->progeny[k], &gparts[offset]);
    }
  }

  /* Return the total number of linked particles. */
  return c->gcount;
}

199
200
201
202
203
204
205
206
207
/**
 * @brief Pack the data of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
 * @param pc Pointer to an array of packed cells in which the
 *      cells will be packed.
 *
 * @return The number of packed cells.
 */
208
209
int cell_pack(struct cell *c, struct pcell *pc) {

210
211
#ifdef WITH_MPI

212
213
  /* Start by packing the data of the current cell. */
  pc->h_max = c->h_max;
214
215
  pc->ti_end_min = c->ti_end_min;
  pc->ti_end_max = c->ti_end_max;
216
  pc->ti_old = c->ti_old;
217
  pc->count = c->count;
218
  pc->gcount = c->gcount;
219
220
221
  c->tag = pc->tag = atomic_inc(&cell_next_tag) % cell_max_tag;

  /* Fill in the progeny, depth-first recursion. */
Pedro Gonnet's avatar
Pedro Gonnet committed
222
223
  int count = 1;
  for (int k = 0; k < 8; k++)
224
225
226
227
228
229
230
    if (c->progeny[k] != NULL) {
      pc->progeny[k] = count;
      count += cell_pack(c->progeny[k], &pc[count]);
    } else
      pc->progeny[k] = -1;

  /* Return the number of packed cells used. */
231
232
  c->pcell_size = count;
  return count;
233
234
235
236
237

#else
  error("SWIFT was not compiled with MPI support.");
  return 0;
#endif
238
239
}

240
241
242
243
244
245
246
247
/**
 * @brief Pack the time information of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
 * @param ti_ends (output) The time information we pack into
 *
 * @return The number of packed cells.
 */
248
int cell_pack_ti_ends(struct cell *c, integertime_t *ti_ends) {
249

250
251
#ifdef WITH_MPI

252
253
  /* Pack this cell's data. */
  ti_ends[0] = c->ti_end_min;
254

255
256
257
258
259
260
261
262
263
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
      count += cell_pack_ti_ends(c->progeny[k], &ti_ends[count]);
    }

  /* Return the number of packed values. */
  return count;
264
265
266
267
268

#else
  error("SWIFT was not compiled with MPI support.");
  return 0;
#endif
269
270
}

271
272
273
274
/**
 * @brief Unpack the time information of a given cell and its sub-cells.
 *
 * @param c The #cell
275
 * @param ti_ends The time information to unpack
276
277
278
 *
 * @return The number of cells created.
 */
279
int cell_unpack_ti_ends(struct cell *c, integertime_t *ti_ends) {
280

281
282
#ifdef WITH_MPI

283
284
  /* Unpack this cell's data. */
  c->ti_end_min = ti_ends[0];
285

286
287
288
289
290
291
292
293
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
      count += cell_unpack_ti_ends(c->progeny[k], &ti_ends[count]);
    }

  /* Return the number of packed values. */
294
  return count;
295
296
297
298
299

#else
  error("SWIFT was not compiled with MPI support.");
  return 0;
#endif
300
}
301

302
/**
303
 * @brief Lock a cell for access to its array of #part and hold its parents.
304
305
 *
 * @param c The #cell.
306
 * @return 0 on success, 1 on failure
307
 */
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
int cell_locktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
  if (c->hold || lock_trylock(&c->lock) != 0) {
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
  if (c->hold) {

    /* Unlock this cell. */
    if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Climb up the tree and lock/hold/unlock. */
Pedro Gonnet's avatar
Pedro Gonnet committed
330
  struct cell *finger;
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

    /* Lock this cell. */
    if (lock_trylock(&finger->lock) != 0) break;

    /* Increment the hold. */
    atomic_inc(&finger->hold);

    /* Unlock the cell. */
    if (lock_unlock(&finger->lock) != 0) error("Failed to unlock cell.");
  }

  /* If we reached the top of the tree, we're done. */
  if (finger == NULL) {
    TIMER_TOC(timer_locktree);
    return 0;
  }

  /* Otherwise, we hit a snag. */
  else {

    /* Undo the holds up to finger. */
Pedro Gonnet's avatar
Pedro Gonnet committed
353
354
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
355
      atomic_dec(&finger2->hold);
356
357
358
359
360
361
362
363
364
365

    /* Unlock this cell. */
    if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }
}

366
367
368
369
370
371
/**
 * @brief Lock a cell for access to its array of #gpart and hold its parents.
 *
 * @param c The #cell.
 * @return 0 on success, 1 on failure
 */
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
int cell_glocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
  if (c->ghold || lock_trylock(&c->glock) != 0) {
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
  if (c->ghold) {

    /* Unlock this cell. */
    if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Climb up the tree and lock/hold/unlock. */
Pedro Gonnet's avatar
Pedro Gonnet committed
394
  struct cell *finger;
395
396
397
398
399
400
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

    /* Lock this cell. */
    if (lock_trylock(&finger->glock) != 0) break;

    /* Increment the hold. */
401
    atomic_inc(&finger->ghold);
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416

    /* Unlock the cell. */
    if (lock_unlock(&finger->glock) != 0) error("Failed to unlock cell.");
  }

  /* If we reached the top of the tree, we're done. */
  if (finger == NULL) {
    TIMER_TOC(timer_locktree);
    return 0;
  }

  /* Otherwise, we hit a snag. */
  else {

    /* Undo the holds up to finger. */
Pedro Gonnet's avatar
Pedro Gonnet committed
417
418
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
419
      atomic_dec(&finger2->ghold);
420
421
422
423
424
425
426
427
428

    /* Unlock this cell. */
    if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }
}
429

430
/**
431
 * @brief Unlock a cell's parents for access to #part array.
432
433
434
 *
 * @param c The #cell.
 */
435
436
437
438
439
440
441
442
void cell_unlocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to unlock this cell. */
  if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
443
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
444
    atomic_dec(&finger->hold);
445
446
447
448

  TIMER_TOC(timer_locktree);
}

449
450
451
452
453
/**
 * @brief Unlock a cell's parents for access to #gpart array.
 *
 * @param c The #cell.
 */
454
455
456
457
458
459
460
461
void cell_gunlocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to unlock this cell. */
  if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
462
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
463
    atomic_dec(&finger->ghold);
464
465
466
467

  TIMER_TOC(timer_locktree);
}

468
469
470
471
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
472
473
 * @param parts_offset Offset of the cell parts array relative to the
 *        space's parts array, i.e. c->parts - s->parts.
474
475
 * @param buff A buffer with at least max(c->count, c->gcount) entries,
 *        used for sorting indices.
Peter W. Draper's avatar
Peter W. Draper committed
476
477
 * @param gbuff A buffer with at least max(c->count, c->gcount) entries,
 *        used for sorting indices for the gparts.
478
 */
479
480
void cell_split(struct cell *c, ptrdiff_t parts_offset, struct cell_buff *buff,
                struct cell_buff *gbuff) {
481

Pedro Gonnet's avatar
Pedro Gonnet committed
482
  const int count = c->count, gcount = c->gcount;
483
484
485
  struct part *parts = c->parts;
  struct xpart *xparts = c->xparts;
  struct gpart *gparts = c->gparts;
486
487
488
489
490
491
  const double pivot[3] = {c->loc[0] + c->width[0] / 2,
                           c->loc[1] + c->width[1] / 2,
                           c->loc[2] + c->width[2] / 2};
  int bucket_count[8] = {0, 0, 0, 0, 0, 0, 0, 0};
  int bucket_offset[9];

492
493
494
#ifdef SWIFT_DEBUG_CHECKS
  /* Check that the buffs are OK. */
  for (int k = 0; k < count; k++) {
495
    if (buff[k].x[0] != parts[k].x[0] || buff[k].x[1] != parts[k].x[1] ||
496
        buff[k].x[2] != parts[k].x[2])
497
498
499
      error("Inconsistent buff contents.");
  }
#endif /* SWIFT_DEBUG_CHECKS */
500
501
502

  /* Fill the buffer with the indices. */
  for (int k = 0; k < count; k++) {
503
504
    const int bid = (buff[k].x[0] > pivot[0]) * 4 +
                    (buff[k].x[1] > pivot[1]) * 2 + (buff[k].x[2] > pivot[2]);
505
    bucket_count[bid]++;
506
    buff[k].ind = bid;
507
  }
508

509
510
511
512
513
  /* Set the buffer offsets. */
  bucket_offset[0] = 0;
  for (int k = 1; k <= 8; k++) {
    bucket_offset[k] = bucket_offset[k - 1] + bucket_count[k - 1];
    bucket_count[k - 1] = 0;
514
515
  }

516
517
518
519
  /* Run through the buckets, and swap particles to their correct spot. */
  for (int bucket = 0; bucket < 8; bucket++) {
    for (int k = bucket_offset[bucket] + bucket_count[bucket];
         k < bucket_offset[bucket + 1]; k++) {
520
      int bid = buff[k].ind;
521
522
523
      if (bid != bucket) {
        struct part part = parts[k];
        struct xpart xpart = xparts[k];
524
        struct cell_buff temp_buff = buff[k];
525
526
        while (bid != bucket) {
          int j = bucket_offset[bid] + bucket_count[bid]++;
527
          while (buff[j].ind == bid) {
528
529
530
            j++;
            bucket_count[bid]++;
          }
Pedro Gonnet's avatar
Pedro Gonnet committed
531
532
          memswap(&parts[j], &part, sizeof(struct part));
          memswap(&xparts[j], &xpart, sizeof(struct xpart));
533
534
          memswap(&buff[j], &temp_buff, sizeof(struct cell_buff));
          bid = temp_buff.ind;
535
536
537
        }
        parts[k] = part;
        xparts[k] = xpart;
538
        buff[k] = temp_buff;
539
      }
540
      bucket_count[bid]++;
541
542
543
544
    }
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
545
  for (int k = 0; k < 8; k++) {
546
547
548
    c->progeny[k]->count = bucket_count[k];
    c->progeny[k]->parts = &c->parts[bucket_offset[k]];
    c->progeny[k]->xparts = &c->xparts[bucket_offset[k]];
549
550
551
  }

  /* Re-link the gparts. */
552
  if (count > 0 && gcount > 0) part_relink_gparts(parts, count, parts_offset);
553

554
#ifdef SWIFT_DEBUG_CHECKS
555
  /* Check that the buffs are OK. */
556
  for (int k = 1; k < count; k++) {
557
558
    if (buff[k].ind < buff[k - 1].ind) error("Buff not sorted.");
    if (buff[k].x[0] != parts[k].x[0] || buff[k].x[1] != parts[k].x[1] ||
559
        buff[k].x[2] != parts[k].x[2])
560
561
562
      error("Inconsistent buff contents (k=%i).", k);
  }

563
  /* Verify that _all_ the parts have been assigned to a cell. */
564
565
566
567
568
569
570
571
  for (int k = 1; k < 8; k++)
    if (&c->progeny[k - 1]->parts[c->progeny[k - 1]->count] !=
        c->progeny[k]->parts)
      error("Particle sorting failed (internal consistency).");
  if (c->progeny[0]->parts != c->parts)
    error("Particle sorting failed (left edge).");
  if (&c->progeny[7]->parts[c->progeny[7]->count] != &c->parts[count])
    error("Particle sorting failed (right edge).");
572
573

  /* Verify a few sub-cells. */
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
  for (int k = 0; k < c->progeny[0]->count; k++)
    if (c->progeny[0]->parts[k].x[0] > pivot[0] ||
        c->progeny[0]->parts[k].x[1] > pivot[1] ||
        c->progeny[0]->parts[k].x[2] > pivot[2])
      error("Sorting failed (progeny=0).");
  for (int k = 0; k < c->progeny[1]->count; k++)
    if (c->progeny[1]->parts[k].x[0] > pivot[0] ||
        c->progeny[1]->parts[k].x[1] > pivot[1] ||
        c->progeny[1]->parts[k].x[2] <= pivot[2])
      error("Sorting failed (progeny=1).");
  for (int k = 0; k < c->progeny[2]->count; k++)
    if (c->progeny[2]->parts[k].x[0] > pivot[0] ||
        c->progeny[2]->parts[k].x[1] <= pivot[1] ||
        c->progeny[2]->parts[k].x[2] > pivot[2])
      error("Sorting failed (progeny=2).");
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
  for (int k = 0; k < c->progeny[3]->count; k++)
    if (c->progeny[3]->parts[k].x[0] > pivot[0] ||
        c->progeny[3]->parts[k].x[1] <= pivot[1] ||
        c->progeny[3]->parts[k].x[2] <= pivot[2])
      error("Sorting failed (progeny=3).");
  for (int k = 0; k < c->progeny[4]->count; k++)
    if (c->progeny[4]->parts[k].x[0] <= pivot[0] ||
        c->progeny[4]->parts[k].x[1] > pivot[1] ||
        c->progeny[4]->parts[k].x[2] > pivot[2])
      error("Sorting failed (progeny=4).");
  for (int k = 0; k < c->progeny[5]->count; k++)
    if (c->progeny[5]->parts[k].x[0] <= pivot[0] ||
        c->progeny[5]->parts[k].x[1] > pivot[1] ||
        c->progeny[5]->parts[k].x[2] <= pivot[2])
      error("Sorting failed (progeny=5).");
  for (int k = 0; k < c->progeny[6]->count; k++)
    if (c->progeny[6]->parts[k].x[0] <= pivot[0] ||
        c->progeny[6]->parts[k].x[1] <= pivot[1] ||
        c->progeny[6]->parts[k].x[2] > pivot[2])
      error("Sorting failed (progeny=6).");
  for (int k = 0; k < c->progeny[7]->count; k++)
    if (c->progeny[7]->parts[k].x[0] <= pivot[0] ||
        c->progeny[7]->parts[k].x[1] <= pivot[1] ||
        c->progeny[7]->parts[k].x[2] <= pivot[2])
      error("Sorting failed (progeny=7).");
614
#endif
615
616

  /* Now do the same song and dance for the gparts. */
617
618
619
620
  for (int k = 0; k < 8; k++) bucket_count[k] = 0;

  /* Fill the buffer with the indices. */
  for (int k = 0; k < gcount; k++) {
621
622
    const int bid = (gbuff[k].x[0] > pivot[0]) * 4 +
                    (gbuff[k].x[1] > pivot[1]) * 2 + (gbuff[k].x[2] > pivot[2]);
623
    bucket_count[bid]++;
624
    gbuff[k].ind = bid;
625
  }
626
627
628
629
630
631

  /* Set the buffer offsets. */
  bucket_offset[0] = 0;
  for (int k = 1; k <= 8; k++) {
    bucket_offset[k] = bucket_offset[k - 1] + bucket_count[k - 1];
    bucket_count[k - 1] = 0;
632
633
  }

634
635
636
637
  /* Run through the buckets, and swap particles to their correct spot. */
  for (int bucket = 0; bucket < 8; bucket++) {
    for (int k = bucket_offset[bucket] + bucket_count[bucket];
         k < bucket_offset[bucket + 1]; k++) {
638
      int bid = gbuff[k].ind;
639
640
      if (bid != bucket) {
        struct gpart gpart = gparts[k];
641
        struct cell_buff temp_buff = gbuff[k];
642
643
        while (bid != bucket) {
          int j = bucket_offset[bid] + bucket_count[bid]++;
644
          while (gbuff[j].ind == bid) {
645
646
647
            j++;
            bucket_count[bid]++;
          }
Pedro Gonnet's avatar
Pedro Gonnet committed
648
          memswap(&gparts[j], &gpart, sizeof(struct gpart));
649
650
          memswap(&gbuff[j], &temp_buff, sizeof(struct cell_buff));
          bid = temp_buff.ind;
651
652
        }
        gparts[k] = gpart;
653
        gbuff[k] = temp_buff;
654
      }
655
      bucket_count[bid]++;
656
657
658
659
    }
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
660
  for (int k = 0; k < 8; k++) {
661
662
    c->progeny[k]->gcount = bucket_count[k];
    c->progeny[k]->gparts = &c->gparts[bucket_offset[k]];
663
664
665
  }

  /* Re-link the parts. */
666
667
  if (count > 0 && gcount > 0)
    part_relink_parts(gparts, gcount, parts - parts_offset);
668
}
669

670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
/**
 * @brief Sanitizes the smoothing length values of cells by setting large
 * outliers to more sensible values.
 *
 * We compute the mean and standard deviation of the smoothing lengths in
 * logarithmic space and limit values to mean + 4 sigma.
 *
 * @param c The cell.
 */
void cell_sanitize(struct cell *c) {

  const int count = c->count;
  struct part *parts = c->parts;

  /* First collect some statistics */
  float h_mean = 0.f, h_mean2 = 0.f;
  float h_min = FLT_MAX, h_max = 0.f;
  for (int i = 0; i < count; ++i) {

689
    const float h = logf(parts[i].h);
690
691
692
693
694
695
696
697
    h_mean += h;
    h_mean2 += h * h;
    h_max = max(h_max, h);
    h_min = min(h_min, h);
  }
  h_mean /= count;
  h_mean2 /= count;
  const float h_var = h_mean2 - h_mean * h_mean;
698
  const float h_std = (h_var > 0.f) ? sqrtf(h_var) : 0.1f * h_mean;
699
700

  /* Choose a cut */
701
  const float h_limit = expf(h_mean + 4.f * h_std);
702
703

  /* Be verbose this is not innocuous */
704
705
  message("Cell properties: h_min= %f h_max= %f geometric mean= %f.",
          expf(h_min), expf(h_max), expf(h_mean));
706
707
708

  if (c->h_max > h_limit) {

709
    message("Smoothing lengths will be limited to (mean + 4sigma)= %f.",
710
711
712
713
714
715
            h_limit);

    /* Apply the cut */
    for (int i = 0; i < count; ++i) parts->h = min(parts[i].h, h_limit);

    c->h_max = h_limit;
716
717
718
719

  } else {

    message("Smoothing lengths will not be limited.");
720
721
722
  }
}

723
/**
724
 * @brief Converts hydro quantities to a valid state after the initial density
725
 * calculation
726
727
728
729
730
731
732
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_convert_hydro(struct cell *c, void *data) {

  struct part *p = c->parts;
733
  struct xpart *xp = c->xparts;
734
735

  for (int i = 0; i < c->count; ++i) {
736
    hydro_convert_quantities(&p[i], &xp[i]);
737
738
739
  }
}

Matthieu Schaller's avatar
Matthieu Schaller committed
740
741
742
743
744
745
/**
 * @brief Cleans the links in a given cell.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
746
void cell_clean_links(struct cell *c, void *data) {
Matthieu Schaller's avatar
Matthieu Schaller committed
747
  c->density = NULL;
748
  c->gradient = NULL;
Matthieu Schaller's avatar
Matthieu Schaller committed
749
  c->force = NULL;
750
  c->grav = NULL;
Matthieu Schaller's avatar
Matthieu Schaller committed
751
}
752

753
754
755
756
757
758
759
760
761
762
/**
 * @brief Checks that a cell is at the current point in time
 *
 * Calls error() if the cell is not at the current time.
 *
 * @param c Cell to act upon
 * @param data The current time on the integer time-line
 */
void cell_check_drift_point(struct cell *c, void *data) {

763
  integertime_t ti_current = *(integertime_t *)data;
764

765
  if (c->ti_old != ti_current && c->nodeID == engine_rank)
766
    error("Cell in an incorrect time-zone! c->ti_old=%lld ti_current=%lld",
767
768
769
          c->ti_old, ti_current);
}

770
771
772
773
774
775
776
777
778
779
780
781
/**
 * @brief Checks whether the cells are direct neighbours ot not. Both cells have
 * to be of the same size
 *
 * @param ci First #cell.
 * @param cj Second #cell.
 *
 * @todo Deal with periodicity.
 */
int cell_are_neighbours(const struct cell *restrict ci,
                        const struct cell *restrict cj) {

Matthieu Schaller's avatar
Matthieu Schaller committed
782
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
783
  if (ci->width[0] != cj->width[0]) error("Cells of different size !");
784
785
786
#endif

  /* Maximum allowed distance */
787
788
  const double min_dist =
      1.2 * ci->width[0]; /* 1.2 accounts for rounding errors */
789
790
791
792
793

  /* (Manhattan) Distance between the cells */
  for (int k = 0; k < 3; k++) {
    const double center_i = ci->loc[k];
    const double center_j = cj->loc[k];
794
    if (fabs(center_i - center_j) > min_dist) return 0;
795
796
797
798
799
  }

  return 1;
}

800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
/**
 * @brief Computes the multi-pole brutally and compare to the
 * recursively computed one.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_check_multipole(struct cell *c, void *data) {

  struct multipole ma;

  if (c->gcount > 0) {

    /* Brute-force calculation */
    multipole_init(&ma, c->gparts, c->gcount);

    /* Compare with recursive one */
    struct multipole mb = c->multipole;

    if (fabsf(ma.mass - mb.mass) / fabsf(ma.mass + mb.mass) > 1e-5)
      error("Multipole masses are different (%12.15e vs. %12.15e)", ma.mass,
            mb.mass);

    for (int k = 0; k < 3; ++k)
824
      if (fabs(ma.CoM[k] - mb.CoM[k]) / fabs(ma.CoM[k] + mb.CoM[k]) > 1e-5)
825
826
827
        error("Multipole CoM are different (%12.15e vs. %12.15e", ma.CoM[k],
              mb.CoM[k]);

828
#if const_gravity_multipole_order >= 2
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
    if (fabsf(ma.I_xx - mb.I_xx) / fabsf(ma.I_xx + mb.I_xx) > 1e-5 &&
        ma.I_xx > 1e-9)
      error("Multipole I_xx are different (%12.15e vs. %12.15e)", ma.I_xx,
            mb.I_xx);
    if (fabsf(ma.I_yy - mb.I_yy) / fabsf(ma.I_yy + mb.I_yy) > 1e-5 &&
        ma.I_yy > 1e-9)
      error("Multipole I_yy are different (%12.15e vs. %12.15e)", ma.I_yy,
            mb.I_yy);
    if (fabsf(ma.I_zz - mb.I_zz) / fabsf(ma.I_zz + mb.I_zz) > 1e-5 &&
        ma.I_zz > 1e-9)
      error("Multipole I_zz are different (%12.15e vs. %12.15e)", ma.I_zz,
            mb.I_zz);
    if (fabsf(ma.I_xy - mb.I_xy) / fabsf(ma.I_xy + mb.I_xy) > 1e-5 &&
        ma.I_xy > 1e-9)
      error("Multipole I_xy are different (%12.15e vs. %12.15e)", ma.I_xy,
            mb.I_xy);
    if (fabsf(ma.I_xz - mb.I_xz) / fabsf(ma.I_xz + mb.I_xz) > 1e-5 &&
        ma.I_xz > 1e-9)
      error("Multipole I_xz are different (%12.15e vs. %12.15e)", ma.I_xz,
            mb.I_xz);
    if (fabsf(ma.I_yz - mb.I_yz) / fabsf(ma.I_yz + mb.I_yz) > 1e-5 &&
        ma.I_yz > 1e-9)
      error("Multipole I_yz are different (%12.15e vs. %12.15e)", ma.I_yz,
            mb.I_yz);
853
#endif
854
  }
855
856
}

857
/**
858
 * @brief Frees up the memory allocated for this #cell.
859
 *
860
 * @param c The #cell.
861
 */
862
863
864
865
866
867
868
void cell_clean(struct cell *c) {

  free(c->sort);

  /* Recurse */
  for (int k = 0; k < 8; k++)
    if (c->progeny[k]) cell_clean(c->progeny[k]);
869
}
870
871
872
873
874

/**
 * @brief Checks whether a given cell needs drifting or not.
 *
 * @param c the #cell.
875
 * @param e The #engine (holding current time information).
876
877
878
 *
 * @return 1 If the cell needs drifting, 0 otherwise.
 */
879
int cell_is_drift_needed(struct cell *c, const struct engine *e) {
880
881

  /* Do we have at least one active particle in the cell ?*/
882
  if (cell_is_active(c, e)) return 1;
883
884
885
886
887
888
889

  /* Loop over the pair tasks that involve this cell */
  for (struct link *l = c->density; l != NULL; l = l->next) {

    if (l->t->type != task_type_pair && l->t->type != task_type_sub_pair)
      continue;

890
891
892
    /* Is the other cell in the pair active ? */
    if ((l->t->ci == c && cell_is_active(l->t->cj, e)) ||
        (l->t->cj == c && cell_is_active(l->t->ci, e)))
893
      return 1;
894
895
896
897
898
  }

  /* No neighbouring cell has active particles. Drift not necessary */
  return 0;
}
899
900
901
902
903
904

/**
 * @brief Un-skips all the tasks associated with a given cell and checks
 * if the space needs to be rebuilt.
 *
 * @param c the #cell.
Peter W. Draper's avatar
Peter W. Draper committed
905
 * @param s the #scheduler.
906
907
908
 *
 * @return 1 If the space needs rebuilding. 0 otherwise.
 */
909
int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
910

911
912
913
914
#ifdef WITH_MPI
  struct engine *e = s->space->e;
#endif

915
  int rebuild = 0;
Matthieu Schaller's avatar
Matthieu Schaller committed
916
  
917
918
919
920
921
  /* Un-skip the density tasks involved with this cell. */
  for (struct link *l = c->density; l != NULL; l = l->next) {
    struct task *t = l->t;
    const struct cell *ci = t->ci;
    const struct cell *cj = t->cj;
922
    scheduler_activate(s, t);
923
924
925
926
927

    /* Set the correct sorting flags */
    if (t->type == task_type_pair) {
      if (!(ci->sorted & (1 << t->flags))) {
        atomic_or(&ci->sorts->flags, (1 << t->flags));
928
        scheduler_activate(s, ci->sorts);
929
930
931
      }
      if (!(cj->sorted & (1 << t->flags))) {
        atomic_or(&cj->sorts->flags, (1 << t->flags));
932
        scheduler_activate(s, cj->sorts);
933
934
935
936
937
938
939
940
941
      }
    }

    /* Check whether there was too much particle motion */
    if (t->type == task_type_pair || t->type == task_type_sub_pair) {
      if (t->tight &&
          (max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin ||
           ci->dx_max > space_maxreldx * ci->h_max ||
           cj->dx_max > space_maxreldx * cj->h_max))
942
        rebuild = 1;
943
944

#ifdef WITH_MPI
945
      /* Activate the send/recv flags. */
946
      if (ci->nodeID != engine_rank) {
947
948

        /* Activate the tasks to recv foreign cell ci's data. */
949
        scheduler_activate(s, ci->recv_xv);
950
951
952
953
        if (cell_is_active(ci, e)) {
          scheduler_activate(s, ci->recv_rho);
          scheduler_activate(s, ci->recv_ti);
        }
954
955
956
957

        /* Look for the local cell cj's send tasks. */
        struct link *l = NULL;
        for (l = cj->send_xv; l != NULL && l->t->cj->nodeID != ci->nodeID;
958
959
             l = l->next)
          ;
960
        if (l == NULL) error("Missing link to send_xv task.");
961
        scheduler_activate(s, l->t);
Matthieu Schaller's avatar
Matthieu Schaller committed
962

Matthieu Schaller's avatar
Matthieu Schaller committed
963
964
965
966
        if (cj->super->drift)
          scheduler_activate(s, cj->super->drift);
        else
          error("Drift task missing !");
967

968
969
970
971
972
973
974
975
976
977
978
979
980
        if (cell_is_active(cj, e)) {
          for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
               l = l->next)
            ;
          if (l == NULL) error("Missing link to send_rho task.");
          scheduler_activate(s, l->t);

          for (l = cj->send_ti; l != NULL && l->t->cj->nodeID != ci->nodeID;
               l = l->next)
            ;
          if (l == NULL) error("Missing link to send_ti task.");
          scheduler_activate(s, l->t);
        }
981

982
      } else if (cj->nodeID != engine_rank) {
983
984

        /* Activate the tasks to recv foreign cell cj's data. */
985
        scheduler_activate(s, cj->recv_xv);
986
987
988
989
        if (cell_is_active(cj, e)) {
          scheduler_activate(s, cj->recv_rho);
          scheduler_activate(s, cj->recv_ti);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
990

991
992
993
        /* Look for the local cell ci's send tasks. */
        struct link *l = NULL;
        for (l = ci->send_xv; l != NULL && l->t->cj->nodeID != cj->nodeID;
994
995
             l = l->next)
          ;
996
        if (l == NULL) error("Missing link to send_xv task.");
997
        scheduler_activate(s, l->t);
Matthieu Schaller's avatar
Matthieu Schaller committed
998

Matthieu Schaller's avatar
Matthieu Schaller committed
999
1000
1001
1002
        if (ci->super->drift)
          scheduler_activate(s, ci->super->drift);
        else
          error("Drift task missing !");
1003

1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
        if (cell_is_active(ci, e)) {
          for (l = ci->send_rho; l != NULL && l->t->cj->nodeID != cj->nodeID;
               l = l->next)
            ;
          if (l == NULL) error("Missing link to send_rho task.");
          scheduler_activate(s, l->t);

          for (l = ci->send_ti; l != NULL && l->t->cj->nodeID != cj->nodeID;
               l = l->next)
            ;
          if (l == NULL) error("Missing link to send_ti task.");
          scheduler_activate(s, l->t);
        }
1017
1018
1019
1020
1021
1022
      }
#endif
    }
  }

  /* Unskip all the other task types. */
1023
  for (struct link *l = c->gradient; l != NULL; l = l->next)
1024
    scheduler_activate(s, l->t);
1025
  for (struct link *l = c->force; l != NULL; l = l->next)
1026
    scheduler_activate(s, l->t);
1027
  for (struct link *l = c->grav; l != NULL; l = l->next)
1028
1029
1030
1031
    scheduler_activate(s, l->t);
  if (c->extra_ghost != NULL) scheduler_activate(s, c->extra_ghost);
  if (c->ghost != NULL) scheduler_activate(s, c->ghost);
  if (c->init != NULL) scheduler_activate(s, c->init);
Matthieu Schaller's avatar
Matthieu Schaller committed
1032
  if (c->drift != NULL) scheduler_activate(s, c->drift);
1033
1034
  if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
  if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
1035
  if (c->timestep != NULL) scheduler_activate(s, c->timestep);
1036
1037
  if (c->cooling != NULL) scheduler_activate(s, c->cooling);
  if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
1038

1039
  return rebuild;
1040
}
1041

1042
1043
1044
1045
1046
1047
1048
1049
1050
/**
 * @brief Set the super-cell pointers for all cells in a hierarchy.
 *
 * @param c The top-level #cell to play with.
 * @param super Pointer to the deepest cell with tasks in this part of the tree.
 */
void cell_set_super(struct cell *c, struct cell *super) {

  /* Are we in a cell with some kind of self/pair task ? */
1051
  if (super == NULL && c->nr_tasks > 0) super = c;
1052
1053
1054
1055

  /* Set the super-cell */
  c->super = super;

1056
1057
  /* Recurse */
  if (c->split)
1058
1059
1060
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) cell_set_super(c->progeny[k], super);
}
1061

1062
1063
1064
1065
1066
1067
/**
 * @brief Recursively drifts all particles and g-particles in a cell hierarchy.
 *
 * @param c The #cell.
 * @param e The #engine (to get ti_current).
 */
1068
void cell_drift(struct cell *c, const struct engine *e) {
1069
1070

  const double timeBase = e->timeBase;
1071
1072
  const integertime_t ti_old = c->ti_old;
  const integertime_t ti_current = e->ti_current;
1073
1074
1075
1076
1077
1078
1079
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;

  /* Drift from the last time the cell was drifted to the current time */
  const double dt = (ti_current - ti_old) * timeBase;
  float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
1080

1081
  /* if (c->gcount > 0) */
Matthieu Schaller's avatar
Matthieu Schaller committed
1082
1083
  /*   message("dt=%e, ti_old=%lld ti_current=%lld", dt, c->ti_old,
   * e->ti_current); */
1084

1085
  /* Check that we are actually going to move forward. */
Matthieu Schaller's avatar
Matthieu Schaller committed
1086
  if (ti_current < ti_old) error("Attempt to drift to the past");
1087

1088
  /* Are we not in a leaf ? */
1089
1090
1091
1092
1093
1094
  if (c->split) {

    /* Loop over the progeny and collect their data. */
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
        struct cell *cp = c->progeny[k];
1095
1096
        cell_drift(cp, e);
        dx_max = max(dx_max, cp->dx_max);
1097
1098
        h_max = max(h_max, cp->h_max);
      }
1099

Matthieu Schaller's avatar
Matthieu Schaller committed
1100
  } else if (ti_current > ti_old) {
Matthieu Schaller's avatar