cell.c 27.7 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 "atomic.h"
51
#include "error.h"
52
#include "gravity.h"
53
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
54
#include "hydro_properties.h"
55
#include "scheduler.h"
56
57
#include "space.h"
#include "timers.h"
58

59
60
61
/* Global variables. */
int cell_next_tag = 0;

62
63
64
65
66
/**
 * @brief Get the size of the cell subtree.
 *
 * @param c The #cell.
 */
67
int cell_getsize(struct cell *c) {
68

Pedro Gonnet's avatar
Pedro Gonnet committed
69
70
  /* Number of cells in this subtree. */
  int count = 1;
71

72
73
  /* Sum up the progeny if split. */
  if (c->split)
Pedro Gonnet's avatar
Pedro Gonnet committed
74
    for (int k = 0; k < 8; k++)
75
76
77
78
79
80
81
      if (c->progeny[k] != NULL) count += cell_getsize(c->progeny[k]);

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

/**
82
83
84
85
86
87
88
89
 * @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.
 */
90
91
int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {

92
93
#ifdef WITH_MPI

94
95
  /* Unpack the current pcell. */
  c->h_max = pc->h_max;
96
97
  c->ti_end_min = pc->ti_end_min;
  c->ti_end_max = pc->ti_end_max;
98
  c->count = pc->count;
99
  c->gcount = pc->gcount;
100
  c->tag = pc->tag;
Matthieu Schaller's avatar
Matthieu Schaller committed
101

102
103
  /* Number of new cells created. */
  int count = 1;
104
105

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
106
  for (int k = 0; k < 8; k++)
107
    if (pc->progeny[k] >= 0) {
Pedro Gonnet's avatar
Pedro Gonnet committed
108
      struct cell *temp = space_getcell(s);
109
      temp->count = 0;
110
      temp->gcount = 0;
111
112
113
      temp->loc[0] = c->loc[0];
      temp->loc[1] = c->loc[1];
      temp->loc[2] = c->loc[2];
114
115
116
      temp->width[0] = c->width[0] / 2;
      temp->width[1] = c->width[1] / 2;
      temp->width[2] = c->width[2] / 2;
117
      temp->dmin = c->dmin / 2;
118
119
120
      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];
121
122
      temp->depth = c->depth + 1;
      temp->split = 0;
123
      temp->dx_max = 0.f;
124
125
126
127
128
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
      count += cell_unpack(&pc[pc->progeny[k]], temp, s);
129
130
    }

131
  /* Return the total number of unpacked cells. */
132
  c->pcell_size = count;
133
  return count;
134
135
136
137
138

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

141
/**
142
 * @brief Link the cells recursively to the given #part array.
143
144
145
146
147
148
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */
149
int cell_link_parts(struct cell *c, struct part *parts) {
150

151
152
153
  c->parts = parts;

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
154
155
156
157
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
158
        offset += cell_link_parts(c->progeny[k], &parts[offset]);
Pedro Gonnet's avatar
Pedro Gonnet committed
159
160
    }
  }
161

162
  /* Return the total number of linked particles. */
163
164
  return c->count;
}
165

166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
/**
 * @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;
}

191
192
193
194
195
196
197
198
199
/**
 * @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.
 */
200
201
int cell_pack(struct cell *c, struct pcell *pc) {

202
203
#ifdef WITH_MPI

204
205
  /* Start by packing the data of the current cell. */
  pc->h_max = c->h_max;
206
207
  pc->ti_end_min = c->ti_end_min;
  pc->ti_end_max = c->ti_end_max;
208
  pc->count = c->count;
209
  pc->gcount = c->gcount;
210
211
212
  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
213
214
  int count = 1;
  for (int k = 0; k < 8; k++)
215
216
217
218
219
220
221
    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. */
222
223
  c->pcell_size = count;
  return count;
224
225
226
227
228

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

231
232
233
234
235
236
237
238
/**
 * @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.
 */
239
240
int cell_pack_ti_ends(struct cell *c, int *ti_ends) {

241
242
#ifdef WITH_MPI

243
244
  /* Pack this cell's data. */
  ti_ends[0] = c->ti_end_min;
245

246
247
248
249
250
251
252
253
254
  /* 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;
255
256
257
258
259

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

262
263
264
265
/**
 * @brief Unpack the time information of a given cell and its sub-cells.
 *
 * @param c The #cell
266
 * @param ti_ends The time information to unpack
267
268
269
 *
 * @return The number of cells created.
 */
270
271
int cell_unpack_ti_ends(struct cell *c, int *ti_ends) {

272
273
#ifdef WITH_MPI

274
275
  /* Unpack this cell's data. */
  c->ti_end_min = ti_ends[0];
276

277
278
279
280
281
282
283
284
  /* 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. */
285
  return count;
286
287
288
289
290

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

293
/**
294
 * @brief Lock a cell for access to its array of #part and hold its parents.
295
296
 *
 * @param c The #cell.
297
 * @return 0 on success, 1 on failure
298
 */
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
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
321
  struct cell *finger;
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
  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
344
345
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
346
      atomic_dec(&finger2->hold);
347
348
349
350
351
352
353
354
355
356

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

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

357
358
359
360
361
362
/**
 * @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
 */
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
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
385
  struct cell *finger;
386
387
388
389
390
391
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

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

    /* Increment the hold. */
392
    atomic_inc(&finger->ghold);
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407

    /* 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
408
409
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
410
      atomic_dec(&finger2->ghold);
411
412
413
414
415
416
417
418
419

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

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

421
/**
422
 * @brief Unlock a cell's parents for access to #part array.
423
424
425
 *
 * @param c The #cell.
 */
426
427
428
429
430
431
432
433
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
434
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
435
    atomic_dec(&finger->hold);
436
437
438
439

  TIMER_TOC(timer_locktree);
}

440
441
442
443
444
/**
 * @brief Unlock a cell's parents for access to #gpart array.
 *
 * @param c The #cell.
 */
445
446
447
448
449
450
451
452
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
453
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
454
    atomic_dec(&finger->ghold);
455
456
457
458

  TIMER_TOC(timer_locktree);
}

459
460
461
462
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
463
464
 * @param parts_offset Offset of the cell parts array relative to the
 *        space's parts array, i.e. c->parts - s->parts.
465
 */
466
void cell_split(struct cell *c, ptrdiff_t parts_offset) {
467

Pedro Gonnet's avatar
Pedro Gonnet committed
468
469
  int i, j;
  const int count = c->count, gcount = c->gcount;
470
471
472
  struct part *parts = c->parts;
  struct xpart *xparts = c->xparts;
  struct gpart *gparts = c->gparts;
473
474
475
476
  int left[8], right[8];
  double pivot[3];

  /* Init the pivots. */
477
  for (int k = 0; k < 3; k++) pivot[k] = c->loc[k] + c->width[k] / 2;
478
479
480
481
482
483
484
485

  /* Split along the x-axis. */
  i = 0;
  j = count - 1;
  while (i <= j) {
    while (i <= count - 1 && parts[i].x[0] <= pivot[0]) i += 1;
    while (j >= 0 && parts[j].x[0] > pivot[0]) j -= 1;
    if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
486
      struct part temp = parts[i];
487
488
      parts[i] = parts[j];
      parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
489
      struct xpart xtemp = xparts[i];
490
491
492
493
      xparts[i] = xparts[j];
      xparts[j] = xtemp;
    }
  }
494
495
496
497
498
499
500
501

#ifdef SWIFT_DEBUG_CHECKS
  for (int k = 0; k <= j; k++)
    if (parts[k].x[0] > pivot[0]) error("cell_split: sorting failed.");
  for (int k = i; k < count; k++)
    if (parts[k].x[0] < pivot[0]) error("cell_split: sorting failed.");
#endif

502
503
504
505
506
507
  left[1] = i;
  right[1] = count - 1;
  left[0] = 0;
  right[0] = j;

  /* Split along the y axis, twice. */
Pedro Gonnet's avatar
Pedro Gonnet committed
508
  for (int k = 1; k >= 0; k--) {
509
510
511
512
513
514
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && parts[i].x[1] <= pivot[1]) i += 1;
      while (j >= left[k] && parts[j].x[1] > pivot[1]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
515
        struct part temp = parts[i];
516
517
        parts[i] = parts[j];
        parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
518
        struct xpart xtemp = xparts[i];
519
520
521
522
        xparts[i] = xparts[j];
        xparts[j] = xtemp;
      }
    }
523
524
525
526
527
528
529
530
531
532
533

#ifdef SWIFT_DEBUG_CHECKS
    for (int kk = left[k]; kk <= j; kk++)
      if (parts[kk].x[1] > pivot[1]) {
        message("ival=[%i,%i], i=%i, j=%i.", left[k], right[k], i, j);
        error("sorting failed (left).");
      }
    for (int kk = i; kk <= right[k]; kk++)
      if (parts[kk].x[1] < pivot[1]) error("sorting failed (right).");
#endif

534
535
536
537
538
539
540
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Split along the z axis, four times. */
Pedro Gonnet's avatar
Pedro Gonnet committed
541
  for (int k = 3; k >= 0; k--) {
542
543
544
545
546
547
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && parts[i].x[2] <= pivot[2]) i += 1;
      while (j >= left[k] && parts[j].x[2] > pivot[2]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
548
        struct part temp = parts[i];
549
550
        parts[i] = parts[j];
        parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
551
        struct xpart xtemp = xparts[i];
552
553
554
555
        xparts[i] = xparts[j];
        xparts[j] = xtemp;
      }
    }
556
557
558
559
560
561
562
563
564
565
566
567
568
569

#ifdef SWIFT_DEBUG_CHECKS
    for (int kk = left[k]; kk <= j; kk++)
      if (parts[kk].x[2] > pivot[2]) {
        message("ival=[%i,%i], i=%i, j=%i.", left[k], right[k], i, j);
        error("sorting failed (left).");
      }
    for (int kk = i; kk <= right[k]; kk++)
      if (parts[kk].x[2] < pivot[2]) {
        message("ival=[%i,%i], i=%i, j=%i.", left[k], right[k], i, j);
        error("sorting failed (right).");
      }
#endif

570
571
572
573
574
575
576
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
577
  for (int k = 0; k < 8; k++) {
578
579
580
581
582
583
    c->progeny[k]->count = right[k] - left[k] + 1;
    c->progeny[k]->parts = &c->parts[left[k]];
    c->progeny[k]->xparts = &c->xparts[left[k]];
  }

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

586
#ifdef SWIFT_DEBUG_CHECKS
587
  /* Verify that _all_ the parts have been assigned to a cell. */
588
589
590
591
592
593
594
595
  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).");
596
597

  /* Verify a few sub-cells. */
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
  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).");
#endif
614
615
616
617
618
619
620
621
622
623

  /* Now do the same song and dance for the gparts. */

  /* Split along the x-axis. */
  i = 0;
  j = gcount - 1;
  while (i <= j) {
    while (i <= gcount - 1 && gparts[i].x[0] <= pivot[0]) i += 1;
    while (j >= 0 && gparts[j].x[0] > pivot[0]) j -= 1;
    if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
624
      struct gpart gtemp = gparts[i];
625
626
      gparts[i] = gparts[j];
      gparts[j] = gtemp;
627
    }
628
629
630
631
632
633
634
  }
  left[1] = i;
  right[1] = gcount - 1;
  left[0] = 0;
  right[0] = j;

  /* Split along the y axis, twice. */
Pedro Gonnet's avatar
Pedro Gonnet committed
635
  for (int k = 1; k >= 0; k--) {
636
637
638
639
640
641
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && gparts[i].x[1] <= pivot[1]) i += 1;
      while (j >= left[k] && gparts[j].x[1] > pivot[1]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
642
        struct gpart gtemp = gparts[i];
643
644
645
646
647
648
649
650
651
652
653
        gparts[i] = gparts[j];
        gparts[j] = gtemp;
      }
    }
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Split along the z axis, four times. */
Pedro Gonnet's avatar
Pedro Gonnet committed
654
  for (int k = 3; k >= 0; k--) {
655
656
657
658
659
660
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && gparts[i].x[2] <= pivot[2]) i += 1;
      while (j >= left[k] && gparts[j].x[2] > pivot[2]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
661
        struct gpart gtemp = gparts[i];
662
663
664
665
666
667
668
669
670
671
672
        gparts[i] = gparts[j];
        gparts[j] = gtemp;
      }
    }
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
673
  for (int k = 0; k < 8; k++) {
674
675
676
677
678
    c->progeny[k]->gcount = right[k] - left[k] + 1;
    c->progeny[k]->gparts = &c->gparts[left[k]];
  }

  /* Re-link the parts. */
679
680
  if (count > 0 && gcount > 0)
    part_relink_parts(gparts, gcount, parts - parts_offset);
681
}
682

683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
/**
 * @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) {

702
    const float h = logf(parts[i].h);
703
704
705
706
707
708
709
710
    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;
711
  const float h_std = (h_var > 0.f) ? sqrtf(h_var) : 0.1f * h_mean;
712
713

  /* Choose a cut */
714
  const float h_limit = expf(h_mean + 4.f * h_std);
715
716

  /* Be verbose this is not innocuous */
717
718
  message("Cell properties: h_min= %f h_max= %f geometric mean= %f.",
          expf(h_min), expf(h_max), expf(h_mean));
719
720
721

  if (c->h_max > h_limit) {

722
    message("Smoothing lengths will be limited to (mean + 4sigma)= %f.",
723
724
725
726
727
728
            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;
729
730
731
732

  } else {

    message("Smoothing lengths will not be limited.");
733
734
735
  }
}

736
/**
737
 * @brief Converts hydro quantities to a valid state after the initial density
738
 * calculation
739
740
741
742
743
744
745
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_convert_hydro(struct cell *c, void *data) {

  struct part *p = c->parts;
746
747

  for (int i = 0; i < c->count; ++i) {
748
749
750
751
    hydro_convert_quantities(&p[i]);
  }
}

Matthieu Schaller's avatar
Matthieu Schaller committed
752
753
754
755
756
757
/**
 * @brief Cleans the links in a given cell.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
758
void cell_clean_links(struct cell *c, void *data) {
Matthieu Schaller's avatar
Matthieu Schaller committed
759
760
  c->density = NULL;
  c->nr_density = 0;
761

762
763
764
  c->gradient = NULL;
  c->nr_gradient = 0;

Matthieu Schaller's avatar
Matthieu Schaller committed
765
766
767
  c->force = NULL;
  c->nr_force = 0;
}
768

769
770
771
772
773
774
775
776
777
778
779
780
/**
 * @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
781
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
782
  if (ci->width[0] != cj->width[0]) error("Cells of different size !");
783
784
785
#endif

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

  /* (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];
    if (fabsf(center_i - center_j) > min_dist) return 0;
  }

  return 1;
}

799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
/**
 * @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)
      if (fabsf(ma.CoM[k] - mb.CoM[k]) / fabsf(ma.CoM[k] + mb.CoM[k]) > 1e-5)
        error("Multipole CoM are different (%12.15e vs. %12.15e", ma.CoM[k],
              mb.CoM[k]);

    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);
  }
852
853
}

854
/**
855
 * @brief Frees up the memory allocated for this #cell.
856
 *
857
 * @param c The #cell.
858
 */
859
860
861
862
863
864
865
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]);
866
}
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886

/**
 * @brief Checks whether a given cell needs drifting or not.
 *
 * @param c the #cell.
 * @param ti_current The current time on the integer time-line.
 *
 * @return 1 If the cell needs drifting, 0 otherwise.
 */
int cell_is_drift_needed(struct cell *c, int ti_current) {

  /* Do we have at least one active particle in the cell ?*/
  if (c->ti_end_min == ti_current) return 1;

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

887
    /* Does the other cell in the pair have an active particle ? */
Matthieu Schaller's avatar
Matthieu Schaller committed
888
889
    if ((l->t->ci == c && l->t->cj->ti_end_min == ti_current) ||
        (l->t->cj == c && l->t->ci->ti_end_min == ti_current))
890
      return 1;
891
892
893
894
895
  }

  /* No neighbouring cell has active particles. Drift not necessary */
  return 0;
}
896
897
898
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.
 *
 * @return 1 If the space needs rebuilding. 0 otherwise.
 */
905
int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
906
907
908
909
910
911

  /* 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;
912
    scheduler_activate(s, t);
913
914
915
916
917

    /* 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));
918
        scheduler_activate(s, ci->sorts);
919
920
921
      }
      if (!(cj->sorted & (1 << t->flags))) {
        atomic_or(&cj->sorts->flags, (1 << t->flags));
922
        scheduler_activate(s, cj->sorts);
923
924
925
926
927
928
929
930
931
932
933
934
      }
    }

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

#ifdef WITH_MPI
935
      /* Activate the send/recv flags. */
936
      if (ci->nodeID != engine_rank) {
937
938

        /* Activate the tasks to recv foreign cell ci's data. */
939
940
941
        scheduler_activate(s, ci->recv_xv);
        scheduler_activate(s, ci->recv_rho);
        scheduler_activate(s, ci->recv_ti);
942
943
944
945

        /* 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;
946
947
             l = l->next)
          ;
948
        if (l == NULL) error("Missing link to send_xv task.");
949
        scheduler_activate(s, l->t);
950
951

        for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
952
953
             l = l->next)
          ;
954
        if (l == NULL) error("Missing link to send_rho task.");
955
        scheduler_activate(s, l->t);
956
957

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

963
      } else if (cj->nodeID != engine_rank) {
964
965

        /* Activate the tasks to recv foreign cell cj's data. */
966
967
968
        scheduler_activate(s, cj->recv_xv);
        scheduler_activate(s, cj->recv_rho);
        scheduler_activate(s, cj->recv_ti);
969
970
971
        /* 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;
972
973
             l = l->next)
          ;
974
        if (l == NULL) error("Missing link to send_xv task.");
975
        scheduler_activate(s, l->t);
976
977

        for (l = ci->send_rho; l != NULL && l->t->cj->nodeID != cj->nodeID;
978
979
             l = l->next)
          ;
980
        if (l == NULL) error("Missing link to send_rho task.");
981
        scheduler_activate(s, l->t);
982
983

        for (l = ci->send_ti; l != NULL && l->t->cj->nodeID != cj->nodeID;
984
985
             l = l->next)
          ;
986
        if (l == NULL) error("Missing link to send_ti task.");
987
        scheduler_activate(s, l->t);
988
989
990
991
992
993
      }
#endif
    }
  }

  /* Unskip all the other task types. */
994
  for (struct link *l = c->gradient; l != NULL; l = l->next)
995
    scheduler_activate(s, l->t);
996
  for (struct link *l = c->force; l != NULL; l = l->next)
997
    scheduler_activate(s, l->t);
998
  for (struct link *l = c->grav; l != NULL; l = l->next)
999
1000
1001
1002
1003
1004
1005
1006
    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);
  if (c->kick != NULL) scheduler_activate(s, c->kick);
  if (c->cooling != NULL) scheduler_activate(s, c->cooling);
  if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
  if (c->grav_external != NULL) scheduler_activate(s, c->grav_external);
1007
1008
1009

  return 0;
}
1010

1011
1012
1013
1014
1015
1016
1017
1018
1019
/**
 * @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 ? */
1020
  if (super == NULL && c->nr_tasks > 0) super = c;
1021
1022
1023
1024

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

1025
1026
  /* Recurse */
  if (c->split)
1027
1028
1029
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) cell_set_super(c->progeny[k], super);
}