cell.c 132 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 "chemistry.h"
53
#include "drift.h"
54
#include "engine.h"
55
#include "error.h"
56
#include "gravity.h"
57
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
58
#include "hydro_properties.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
59
#include "memswap.h"
60
#include "minmax.h"
61
#include "scheduler.h"
62
#include "space.h"
63
#include "space_getsid.h"
Loic Hausammann's avatar
Loic Hausammann committed
64
#include "stars.h"
65
#include "timers.h"
66
#include "tools.h"
67

68
69
70
/* Global variables. */
int cell_next_tag = 0;

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

Pedro Gonnet's avatar
Pedro Gonnet committed
78
79
  /* Number of cells in this subtree. */
  int count = 1;
80

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

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

90
/**
91
 * @brief Link the cells recursively to the given #part array.
92
93
94
95
96
97
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */
98
int cell_link_parts(struct cell *c, struct part *parts) {
99

100
  c->hydro.parts = parts;
101
102

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
103
104
105
106
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
107
        offset += cell_link_parts(c->progeny[k], &parts[offset]);
Pedro Gonnet's avatar
Pedro Gonnet committed
108
109
    }
  }
110

111
  /* Return the total number of linked particles. */
112
  return c->hydro.count;
113
}
114

115
116
117
118
119
120
121
122
123
124
/**
 * @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) {

125
  c->grav.parts = gparts;
126
127
128
129
130
131
132
133
134
135
136

  /* 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. */
137
  return c->grav.count;
138
139
}

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

150
  c->stars.parts = sparts;
151
152
153
154
155
156
157
158
159
160
161

  /* 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_sparts(c->progeny[k], &sparts[offset]);
    }
  }

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

165
166
167
168
169
170
/**
 * @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.
171
172
 * @param with_gravity Are we running with gravity and hence need
 *      to exchange multipoles?
173
174
175
 *
 * @return The number of packed cells.
 */
176
int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
Matthieu Schaller's avatar
Matthieu Schaller committed
177
              const int with_gravity) {
178

179
180
#ifdef WITH_MPI

181
  /* Start by packing the data of the current cell. */
182
183
184
185
186
  pc->hydro.h_max = c->hydro.h_max;
  pc->hydro.ti_end_min = c->hydro.ti_end_min;
  pc->hydro.ti_end_max = c->hydro.ti_end_max;
  pc->grav.ti_end_min = c->grav.ti_end_min;
  pc->grav.ti_end_max = c->grav.ti_end_max;
187
188
  pc->hydro.ti_old_part = c->hydro.ti_old_part;
  pc->grav.ti_old_part = c->grav.ti_old_part;
189
190
  pc->grav.ti_old_multipole = c->grav.ti_old_multipole;
  pc->hydro.count = c->hydro.count;
191
192
  pc->grav.count = c->grav.count;
  pc->stars.count = c->stars.count;
193

194
  /* Copy the Multipole related information */
Matthieu Schaller's avatar
Matthieu Schaller committed
195
  if (with_gravity) {
196
    const struct gravity_tensors *mp = c->grav.multipole;
197

198
199
200
201
202
203
204
205
206
    pc->grav.m_pole = mp->m_pole;
    pc->grav.CoM[0] = mp->CoM[0];
    pc->grav.CoM[1] = mp->CoM[1];
    pc->grav.CoM[2] = mp->CoM[2];
    pc->grav.CoM_rebuild[0] = mp->CoM_rebuild[0];
    pc->grav.CoM_rebuild[1] = mp->CoM_rebuild[1];
    pc->grav.CoM_rebuild[2] = mp->CoM_rebuild[2];
    pc->grav.r_max = mp->r_max;
    pc->grav.r_max_rebuild = mp->r_max_rebuild;
207
208
  }

209
210
211
#ifdef SWIFT_DEBUG_CHECKS
  pc->cellID = c->cellID;
#endif
212
213

  /* Fill in the progeny, depth-first recursion. */
Pedro Gonnet's avatar
Pedro Gonnet committed
214
215
  int count = 1;
  for (int k = 0; k < 8; k++)
216
217
    if (c->progeny[k] != NULL) {
      pc->progeny[k] = count;
218
      count += cell_pack(c->progeny[k], &pc[count], with_gravity);
219
    } else {
220
      pc->progeny[k] = -1;
221
    }
222
223

  /* Return the number of packed cells used. */
224
  c->mpi.pcell_size = count;
225
  return count;
226
227
228
229
230

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

233
234
235
236
237
238
239
240
241
242
243
244
245
/**
 * @brief Pack the tag of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
 * @param tags Pointer to an array of packed tags.
 *
 * @return The number of packed tags.
 */
int cell_pack_tags(const struct cell *c, int *tags) {

#ifdef WITH_MPI

  /* Start by packing the data of the current cell. */
246
  tags[0] = c->mpi.tag;
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_tags(c->progeny[k], &tags[count]);

#ifdef SWIFT_DEBUG_CHECKS
255
  if (c->mpi.pcell_size != count) error("Inconsistent tag and pcell count!");
256
257
258
259
260
261
262
263
264
265
266
#endif  // SWIFT_DEBUG_CHECKS

  /* Return the number of packed tags used. */
  return count;

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

267
268
269
270
271
272
/**
 * @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.
273
274
 * @param with_gravity Are we running with gravity and hence need
 *      to exchange multipoles?
275
276
277
 *
 * @return The number of cells created.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
278
int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
279
                struct space *restrict s, const int with_gravity) {
280
281
282
283

#ifdef WITH_MPI

  /* Unpack the current pcell. */
284
285
286
287
288
  c->hydro.h_max = pc->hydro.h_max;
  c->hydro.ti_end_min = pc->hydro.ti_end_min;
  c->hydro.ti_end_max = pc->hydro.ti_end_max;
  c->grav.ti_end_min = pc->grav.ti_end_min;
  c->grav.ti_end_max = pc->grav.ti_end_max;
289
290
  c->hydro.ti_old_part = pc->hydro.ti_old_part;
  c->grav.ti_old_part = pc->grav.ti_old_part;
291
292
  c->grav.ti_old_multipole = pc->grav.ti_old_multipole;
  c->hydro.count = pc->hydro.count;
293
294
  c->grav.count = pc->grav.count;
  c->stars.count = pc->stars.count;
295
296
297
#ifdef SWIFT_DEBUG_CHECKS
  c->cellID = pc->cellID;
#endif
298

299
  /* Copy the Multipole related information */
Matthieu Schaller's avatar
Matthieu Schaller committed
300
  if (with_gravity) {
301

302
    struct gravity_tensors *mp = c->grav.multipole;
303

304
305
306
307
308
309
310
311
312
    mp->m_pole = pc->grav.m_pole;
    mp->CoM[0] = pc->grav.CoM[0];
    mp->CoM[1] = pc->grav.CoM[1];
    mp->CoM[2] = pc->grav.CoM[2];
    mp->CoM_rebuild[0] = pc->grav.CoM_rebuild[0];
    mp->CoM_rebuild[1] = pc->grav.CoM_rebuild[1];
    mp->CoM_rebuild[2] = pc->grav.CoM_rebuild[2];
    mp->r_max = pc->grav.r_max;
    mp->r_max_rebuild = pc->grav.r_max_rebuild;
313
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
314

315
316
317
318
  /* Number of new cells created. */
  int count = 1;

  /* Fill the progeny recursively, depth-first. */
319
  c->split = 0;
320
321
322
323
  for (int k = 0; k < 8; k++)
    if (pc->progeny[k] >= 0) {
      struct cell *temp;
      space_getcells(s, 1, &temp);
324
      temp->hydro.count = 0;
325
326
      temp->grav.count = 0;
      temp->stars.count = 0;
327
328
329
330
331
332
333
334
335
336
337
338
      temp->loc[0] = c->loc[0];
      temp->loc[1] = c->loc[1];
      temp->loc[2] = c->loc[2];
      temp->width[0] = c->width[0] / 2;
      temp->width[1] = c->width[1] / 2;
      temp->width[2] = c->width[2] / 2;
      temp->dmin = c->dmin / 2;
      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];
      temp->depth = c->depth + 1;
      temp->split = 0;
339
      temp->hydro.dx_max_part = 0.f;
340
      temp->hydro.dx_max_sort = 0.f;
Loic Hausammann's avatar
Loic Hausammann committed
341
      temp->stars.dx_max_part = 0.f;
342
343
344
345
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
346
      count += cell_unpack(&pc[pc->progeny[k]], temp, s, with_gravity);
347
348
349
    }

  /* Return the total number of unpacked cells. */
350
  c->mpi.pcell_size = count;
351
352
353
354
355
356
357
358
  return count;

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

359
360
361
362
363
364
365
366
367
368
369
370
371
/**
 * @brief Unpack the tags of a given cell and its sub-cells.
 *
 * @param tags An array of tags.
 * @param c The #cell in which to unpack the tags.
 *
 * @return The number of tags created.
 */
int cell_unpack_tags(const int *tags, struct cell *restrict c) {

#ifdef WITH_MPI

  /* Unpack the current pcell. */
372
  c->mpi.tag = tags[0];
373
374
375
376
377
378
379
380
381
382
383

  /* Number of new cells created. */
  int count = 1;

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

#ifdef SWIFT_DEBUG_CHECKS
384
  if (c->mpi.pcell_size != count) error("Inconsistent tag and pcell count!");
385
386
387
388
389
390
391
392
393
394
395
#endif  // SWIFT_DEBUG_CHECKS

  /* Return the total number of unpacked tags. */
  return count;

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

396
397
398
399
/**
 * @brief Pack the time information of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
400
 * @param pcells (output) The end-of-timestep information we pack into
401
402
403
 *
 * @return The number of packed cells.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
404
405
int cell_pack_end_step(struct cell *restrict c,
                       struct pcell_step *restrict pcells) {
406

407
408
#ifdef WITH_MPI

409
  /* Pack this cell's data. */
410
411
412
413
  pcells[0].hydro.ti_end_min = c->hydro.ti_end_min;
  pcells[0].hydro.ti_end_max = c->hydro.ti_end_max;
  pcells[0].grav.ti_end_min = c->grav.ti_end_min;
  pcells[0].grav.ti_end_max = c->grav.ti_end_max;
414
  pcells[0].hydro.dx_max_part = c->hydro.dx_max_part;
Loic Hausammann's avatar
Loic Hausammann committed
415
  pcells[0].stars.dx_max_part = c->stars.dx_max_part;
416

417
418
419
420
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
421
      count += cell_pack_end_step(c->progeny[k], &pcells[count]);
422
423
424
425
    }

  /* Return the number of packed values. */
  return count;
426
427
428
429
430

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

433
434
435
436
/**
 * @brief Unpack the time information of a given cell and its sub-cells.
 *
 * @param c The #cell
437
 * @param pcells The end-of-timestep information to unpack
438
439
440
 *
 * @return The number of cells created.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
441
442
int cell_unpack_end_step(struct cell *restrict c,
                         struct pcell_step *restrict pcells) {
443

444
445
#ifdef WITH_MPI

446
  /* Unpack this cell's data. */
447
448
449
450
  c->hydro.ti_end_min = pcells[0].hydro.ti_end_min;
  c->hydro.ti_end_max = pcells[0].hydro.ti_end_max;
  c->grav.ti_end_min = pcells[0].grav.ti_end_min;
  c->grav.ti_end_max = pcells[0].grav.ti_end_max;
451
  c->hydro.dx_max_part = pcells[0].hydro.dx_max_part;
Loic Hausammann's avatar
Loic Hausammann committed
452
  c->stars.dx_max_part = pcells[0].stars.dx_max_part;
453

454
455
456
457
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
458
      count += cell_unpack_end_step(c->progeny[k], &pcells[count]);
459
460
461
    }

  /* Return the number of packed values. */
462
  return count;
463
464
465
466
467

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

470
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
471
472
 * @brief Pack the multipole information of the given cell and all it's
 * sub-cells.
473
474
475
476
477
478
479
 *
 * @param c The #cell.
 * @param pcells (output) The multipole information we pack into
 *
 * @return The number of packed cells.
 */
int cell_pack_multipoles(struct cell *restrict c,
Matthieu Schaller's avatar
Matthieu Schaller committed
480
                         struct gravity_tensors *restrict pcells) {
481
482
483
484

#ifdef WITH_MPI

  /* Pack this cell's data. */
485
  pcells[0] = *c->grav.multipole;
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511

  /* 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_multipoles(c->progeny[k], &pcells[count]);
    }

  /* Return the number of packed values. */
  return count;

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

/**
 * @brief Unpack the multipole information of a given cell and its sub-cells.
 *
 * @param c The #cell
 * @param pcells The multipole information to unpack
 *
 * @return The number of cells created.
 */
int cell_unpack_multipoles(struct cell *restrict c,
Matthieu Schaller's avatar
Matthieu Schaller committed
512
                           struct gravity_tensors *restrict pcells) {
513
514
515
516

#ifdef WITH_MPI

  /* Unpack this cell's data. */
517
  *c->grav.multipole = pcells[0];
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534

  /* 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_multipoles(c->progeny[k], &pcells[count]);
    }

  /* Return the number of packed values. */
  return count;

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

535
/**
536
 * @brief Lock a cell for access to its array of #part and hold its parents.
537
538
 *
 * @param c The #cell.
539
 * @return 0 on success, 1 on failure
540
 */
541
542
543
544
545
int cell_locktree(struct cell *c) {

  TIMER_TIC

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

  /* Did somebody hold this cell in the meantime? */
552
  if (c->hydro.hold) {
553
554

    /* Unlock this cell. */
555
    if (lock_unlock(&c->hydro.lock) != 0) error("Failed to unlock cell.");
556
557
558
559
560
561
562

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

  /* Climb up the tree and lock/hold/unlock. */
Pedro Gonnet's avatar
Pedro Gonnet committed
563
  struct cell *finger;
564
565
566
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

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

    /* Increment the hold. */
570
    atomic_inc(&finger->hydro.hold);
571
572

    /* Unlock the cell. */
573
    if (lock_unlock(&finger->hydro.lock) != 0) error("Failed to unlock cell.");
574
575
576
577
578
579
580
581
582
583
584
585
  }

  /* 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
586
587
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
588
      atomic_dec(&finger2->hydro.hold);
589
590

    /* Unlock this cell. */
591
    if (lock_unlock(&c->hydro.lock) != 0) error("Failed to unlock cell.");
592
593
594
595
596
597
598

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

599
600
601
602
603
604
/**
 * @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
 */
605
606
607
608
609
int cell_glocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
610
  if (c->grav.phold || lock_trylock(&c->grav.plock) != 0) {
611
612
613
614
615
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
616
  if (c->grav.phold) {
617
618

    /* Unlock this cell. */
619
    if (lock_unlock(&c->grav.plock) != 0) error("Failed to unlock cell.");
620
621
622
623
624
625
626

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

  /* Climb up the tree and lock/hold/unlock. */
Pedro Gonnet's avatar
Pedro Gonnet committed
627
  struct cell *finger;
628
629
630
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

    /* Lock this cell. */
631
    if (lock_trylock(&finger->grav.plock) != 0) break;
632
633

    /* Increment the hold. */
634
    atomic_inc(&finger->grav.phold);
635
636

    /* Unlock the cell. */
637
    if (lock_unlock(&finger->grav.plock) != 0) error("Failed to unlock cell.");
638
639
640
641
642
643
644
645
646
647
648
649
  }

  /* 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
650
651
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
652
      atomic_dec(&finger2->grav.phold);
653
654

    /* Unlock this cell. */
655
    if (lock_unlock(&c->grav.plock) != 0) error("Failed to unlock cell.");
656
657
658
659
660
661

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

663
664
665
666
667
668
669
670
671
672
673
/**
 * @brief Lock a cell for access to its #multipole and hold its parents.
 *
 * @param c The #cell.
 * @return 0 on success, 1 on failure
 */
int cell_mlocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
674
  if (c->grav.mhold || lock_trylock(&c->grav.mlock) != 0) {
675
676
677
678
679
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
680
  if (c->grav.mhold) {
681
682

    /* Unlock this cell. */
683
    if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
684
685
686
687
688
689
690
691
692
693
694

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

  /* Climb up the tree and lock/hold/unlock. */
  struct cell *finger;
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

    /* Lock this cell. */
695
    if (lock_trylock(&finger->grav.mlock) != 0) break;
696
697

    /* Increment the hold. */
698
    atomic_inc(&finger->grav.mhold);
699
700

    /* Unlock the cell. */
701
    if (lock_unlock(&finger->grav.mlock) != 0) error("Failed to unlock cell.");
702
703
704
705
706
707
708
709
710
711
712
713
714
715
  }

  /* 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. */
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
716
      atomic_dec(&finger2->grav.mhold);
717
718

    /* Unlock this cell. */
719
    if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
720
721
722
723
724
725
726

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

727
728
729
730
731
732
733
734
735
736
737
/**
 * @brief Lock a cell for access to its array of #spart and hold its parents.
 *
 * @param c The #cell.
 * @return 0 on success, 1 on failure
 */
int cell_slocktree(struct cell *c) {

  TIMER_TIC

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

  /* Did somebody hold this cell in the meantime? */
744
  if (c->stars.hold) {
745
746

    /* Unlock this cell. */
747
    if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
748
749
750
751
752
753
754
755
756
757
758

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

  /* Climb up the tree and lock/hold/unlock. */
  struct cell *finger;
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

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

    /* Increment the hold. */
762
    atomic_inc(&finger->stars.hold);
763
764

    /* Unlock the cell. */
765
    if (lock_unlock(&finger->stars.lock) != 0) error("Failed to unlock cell.");
766
767
768
769
770
771
772
773
774
775
776
777
778
779
  }

  /* 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. */
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
780
      atomic_dec(&finger2->stars.hold);
781
782

    /* Unlock this cell. */
783
    if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
784
785
786
787
788
789
790

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

791
/**
792
 * @brief Unlock a cell's parents for access to #part array.
793
794
795
 *
 * @param c The #cell.
 */
796
797
798
799
800
void cell_unlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
804
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
805
    atomic_dec(&finger->hydro.hold);
806
807
808
809

  TIMER_TOC(timer_locktree);
}

810
811
812
813
814
/**
 * @brief Unlock a cell's parents for access to #gpart array.
 *
 * @param c The #cell.
 */
815
816
817
818
819
void cell_gunlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
823
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
824
    atomic_dec(&finger->grav.phold);
825
826
827
828

  TIMER_TOC(timer_locktree);
}

829
830
831
832
833
834
835
836
837
838
/**
 * @brief Unlock a cell's parents for access to its #multipole.
 *
 * @param c The #cell.
 */
void cell_munlocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to unlock this cell. */
839
  if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
840
841
842

  /* Climb up the tree and unhold the parents. */
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
843
    atomic_dec(&finger->grav.mhold);
844
845
846
847

  TIMER_TOC(timer_locktree);
}

848
849
850
851
852
853
854
855
856
857
/**
 * @brief Unlock a cell's parents for access to #spart array.
 *
 * @param c The #cell.
 */
void cell_sunlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
862
    atomic_dec(&finger->stars.hold);
863
864
865
866

  TIMER_TOC(timer_locktree);
}

867
868
869
870
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
871
 * @param parts_offset Offset of the cell parts array relative to the
872
 *        space's parts array, i.e. c->hydro.parts - s->parts.
873
 * @param sparts_offset Offset of the cell sparts array relative to the
874
875
 *        space's sparts array, i.e. c->stars.parts - s->stars.parts.
 * @param buff A buffer with at least max(c->hydro.count, c->grav.count)
876
 * entries, used for sorting indices.
877
878
879
 * @param sbuff A buffer with at least max(c->stars.count, c->grav.count)
 * entries, used for sorting indices for the sparts.
 * @param gbuff A buffer with at least max(c->hydro.count, c->grav.count)
880
 * entries, used for sorting indices for the gparts.
881
 */
882
883
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
                struct cell_buff *buff, struct cell_buff *sbuff,
884
                struct cell_buff *gbuff) {
885

886
887
  const int count = c->hydro.count, gcount = c->grav.count,
            scount = c->stars.count;
888
889
  struct part *parts = c->hydro.parts;
  struct xpart *xparts = c->hydro.xparts;
890
891
  struct gpart *gparts = c->grav.parts;
  struct spart *sparts = c->stars.parts;
892
893
894
895
896
897
  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];

898
899
900
#ifdef SWIFT_DEBUG_CHECKS
  /* Check that the buffs are OK. */
  for (int k = 0; k < count; k++) {
901
    if (buff[k].x[0] != parts[k].x[0] || buff[k].x[1] != parts[k].x[1] ||
902
        buff[k].x[2] != parts[k].x[2])
903
904
      error("Inconsistent buff contents.");
  }
905
906
907
908
909
910
911
912
913
914
  for (int k = 0; k < gcount; k++) {
    if (gbuff[k].x[0] != gparts[k].x[0] || gbuff[k].x[1] != gparts[k].x[1] ||
        gbuff[k].x[2] != gparts[k].x[2])
      error("Inconsistent gbuff contents.");
  }
  for (int k = 0; k < scount; k++) {
    if (sbuff[k].x[0] != sparts[k].x[0] || sbuff[k].x[1] != sparts[k].x[1] ||
        sbuff[k].x[2] != sparts[k].x[2])
      error("Inconsistent sbuff contents.");
  }
915
#endif /* SWIFT_DEBUG_CHECKS */
916
917
918

  /* Fill the buffer with the indices. */
  for (int k = 0; k < count; k++) {
919
920
    const int bid = (buff[k].x[0] >= pivot[0]) * 4 +
                    (buff[k].x[1] >= pivot[1]) * 2 + (buff[k].x[2] >= pivot[2]);
921
    bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
922
    buff[k].ind = bid;
923
  }
924

925
926
927
928
929
  /* 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;
930
931
  }

932
933
934
935
  /* 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++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
936
      int bid = buff[k].ind;
937
938
939
      if (bid != bucket) {
        struct part part = parts[k];
        struct xpart xpart = xparts[k];
940
        struct cell_buff temp_buff = buff[k];
941
942
        while (bid != bucket) {
          int j = bucket_offset[bid] + bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
943
          while (buff[j].ind == bid) {
944
945
946
            j++;
            bucket_count[bid]++;
          }
Pedro Gonnet's avatar
Pedro Gonnet committed
947
948
          memswap(&parts[j], &part, sizeof(struct part));
          memswap(&xparts[j], &xpart, sizeof(struct xpart));
949
          memswap(&buff[j], &temp_buff, sizeof(struct cell_buff));
950
951
          if (parts[j].gpart)
            parts[j].gpart->id_or_neg_offset = -(j + parts_offset);
952
          bid = temp_buff.ind;
953
954
955
        }
        parts[k] = part;
        xparts[k] = xpart;
956
        buff[k] = temp_buff;
957
958
        if (parts[k].gpart)
          parts[k].gpart->id_or_neg_offset = -(k + parts_offset);
959
      }
960
      bucket_count[bid]++;
961
962
963
964
    }
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
965
  for (int k = 0; k < 8; k++) {
966
967
968
    c->progeny[k]->hydro.count = bucket_count[k];
    c->progeny[k]->hydro.parts = &c->hydro.parts[bucket_offset[k]];
    c->progeny[k]->hydro.xparts = &c->hydro.xparts[bucket_offset[k]];
969
970
  }

971
#ifdef SWIFT_DEBUG_CHECKS
972
973
974
975
976
977
978
979
  /* Check that the buffs are OK. */
  for (int k = 1; k < count; k++) {
    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] ||
        buff[k].x[2] != parts[k].x[2])
      error("Inconsistent buff contents (k=%i).", k);
  }

980
  /* Verify that _all_ the parts have been assigned to a cell. */
981
  for (int k = 1; k < 8; k++)
982
983
    if (&c->progeny[k - 1]->hydro.parts[c->progeny[k - 1]->hydro.count] !=
        c->progeny[k]->hydro.parts)
984
      error("Particle sorting failed (internal consistency).");
985
  if (c->progeny[0]->hydro.parts != c->hydro.parts)
986
    error("Particle sorting failed (left edge).");
987
988
  if (&c->progeny[7]->hydro.parts[c->progeny[7]->hydro.count] !=
      &c->hydro.parts[count])
989
    error("Particle sorting failed (right edge).");
990
991

  /* Verify a few sub-cells. */
992
993
994
995
  for (int k = 0; k < c->progeny[0]->hydro.count; k++)
    if (c->progeny[0]->hydro.parts[k].x[0] >= pivot[0] ||
        c->progeny[0]->hydro.parts[k].x[1] >= pivot[1] ||
        c->progeny[0]->hydro.parts[k].x[2] >= pivot[2])
996
      error("Sorting failed (progeny=0).");
997
998
999
1000
  for (int k = 0; k < c->progeny[1]->hydro.count; k++)
    if (c->progeny[1]->hydro.parts[k].x[0] >= pivot[0] ||
        c->progeny[1]->hydro.parts[k].x[1] >= pivot[1] ||
        c->progeny[1]->hydro.parts[k].x[2] < pivot[2])