cell.c 155 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
#include "tracers.h"
68

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

72
73
74
75
76
int depth_0_counter = 0;
int depth_1_counter = 0;
int depth_2_counter = 0;
int depth_3_counter = 0;

77
78
79
80
81
/**
 * @brief Get the size of the cell subtree.
 *
 * @param c The #cell.
 */
82
int cell_getsize(struct cell *c) {
83

Pedro Gonnet's avatar
Pedro Gonnet committed
84
85
  /* Number of cells in this subtree. */
  int count = 1;
86

87
88
  /* Sum up the progeny if split. */
  if (c->split)
Pedro Gonnet's avatar
Pedro Gonnet committed
89
    for (int k = 0; k < 8; k++)
90
91
92
93
94
95
      if (c->progeny[k] != NULL) count += cell_getsize(c->progeny[k]);

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

96
/**
97
 * @brief Link the cells recursively to the given #part array.
98
99
100
101
102
103
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */
104
105
106
107
108
109
int cell_link_parts(struct cell *c, struct part *parts, int link) {

#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");
#endif
110

111
112
113
114
115
116
117
118
119
  /* Do we have a hydro task at this level? */
  if (c->hydro.density != NULL) {
    link = 1;
  }

  /* Ok, link the particles at this level */
  if (link) {
    c->hydro.parts = parts;
  }
120
121

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
122
  if (c->split) {
123
124
    int count = 0;

Pedro Gonnet's avatar
Pedro Gonnet committed
125
    for (int k = 0; k < 8; k++) {
126
127
128
      if (c->progeny[k] != NULL) {
        count += cell_link_parts(c->progeny[k], &parts[count], link);
      }
Pedro Gonnet's avatar
Pedro Gonnet committed
129
    }
130
131
132
133
134
135
136

#ifdef SWIFT_DEBUG_CHECKS
    if (link && (count != c->hydro.count))
      error("Something is wrong with the foreign part counts.");
#endif

    return count;
Pedro Gonnet's avatar
Pedro Gonnet committed
137
  }
138

139
  /* Return the total number of linked particles. */
140
141
142
143
  if (link)
    return c->hydro.count;
  else
    return 0;
144
}
145

146
147
148
149
150
151
152
153
154
155
/**
 * @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) {

156
  c->grav.parts = gparts;
157
158
159
160
161
162
163
164
165
166
167

  /* 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. */
168
  return c->grav.count;
169
170
}

171
172
173
174
175
176
177
178
179
180
/**
 * @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) {

181
  c->stars.parts = sparts;
182
183
184
185
186
187
188
189
190
191
192

  /* 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. */
193
  return c->stars.count;
194
195
}

196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
int cell_count_parts_for_tasks(const struct cell *c) {

#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Counting foreign particles in a local cell!");
#endif

  /* Do we have a hydro task at this level? */
  if (c->hydro.density != NULL) {
    if (c->depth == 0) ++depth_0_counter;
    if (c->depth == 1) ++depth_1_counter;
    if (c->depth == 2) ++depth_2_counter;
    if (c->depth == 3) ++depth_3_counter;
    return c->hydro.count;
  }

  if (c->split) {
    int count = 0;
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
        count += cell_count_parts_for_tasks(c->progeny[k]);
      }
    }
    return count;
  } else {
    return 0;
  }
}

225
226
227
228
229
230
/**
 * @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.
231
232
 * @param with_gravity Are we running with gravity and hence need
 *      to exchange multipoles?
233
234
235
 *
 * @return The number of packed cells.
 */
236
int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
Matthieu Schaller's avatar
Matthieu Schaller committed
237
              const int with_gravity) {
238

239
240
#ifdef WITH_MPI

241
  /* Start by packing the data of the current cell. */
242
243
244
245
246
  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;
247
  pc->stars.ti_end_min = c->stars.ti_end_min;
248
249
  pc->hydro.ti_old_part = c->hydro.ti_old_part;
  pc->grav.ti_old_part = c->grav.ti_old_part;
250
251
  pc->grav.ti_old_multipole = c->grav.ti_old_multipole;
  pc->hydro.count = c->hydro.count;
252
253
  pc->grav.count = c->grav.count;
  pc->stars.count = c->stars.count;
254
  pc->maxdepth = c->maxdepth;
255

256
  /* Copy the Multipole related information */
Matthieu Schaller's avatar
Matthieu Schaller committed
257
  if (with_gravity) {
258
    const struct gravity_tensors *mp = c->grav.multipole;
259

260
261
262
263
264
265
266
267
268
    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;
269
270
  }

271
272
273
#ifdef SWIFT_DEBUG_CHECKS
  pc->cellID = c->cellID;
#endif
274
275

  /* Fill in the progeny, depth-first recursion. */
Pedro Gonnet's avatar
Pedro Gonnet committed
276
277
  int count = 1;
  for (int k = 0; k < 8; k++)
278
279
    if (c->progeny[k] != NULL) {
      pc->progeny[k] = count;
280
      count += cell_pack(c->progeny[k], &pc[count], with_gravity);
281
    } else {
282
      pc->progeny[k] = -1;
283
    }
284
285

  /* Return the number of packed cells used. */
286
  c->mpi.pcell_size = count;
287
  return count;
288
289
290
291
292

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

295
296
297
298
299
300
301
302
303
304
305
306
307
/**
 * @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. */
308
  tags[0] = c->mpi.tag;
309
310
311
312
313
314
315
316

  /* 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
317
  if (c->mpi.pcell_size != count) error("Inconsistent tag and pcell count!");
318
319
320
321
322
323
324
325
326
327
328
#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
}

329
330
331
332
333
334
/**
 * @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.
335
336
 * @param with_gravity Are we running with gravity and hence need
 *      to exchange multipoles?
337
338
339
 *
 * @return The number of cells created.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
340
int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
341
                struct space *restrict s, const int with_gravity) {
342
343
344
345

#ifdef WITH_MPI

  /* Unpack the current pcell. */
346
347
348
349
350
  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;
351
  c->stars.ti_end_min = pc->stars.ti_end_min;
352
353
  c->hydro.ti_old_part = pc->hydro.ti_old_part;
  c->grav.ti_old_part = pc->grav.ti_old_part;
354
355
  c->grav.ti_old_multipole = pc->grav.ti_old_multipole;
  c->hydro.count = pc->hydro.count;
356
357
  c->grav.count = pc->grav.count;
  c->stars.count = pc->stars.count;
358
359
  c->maxdepth = pc->maxdepth;

360
361
362
#ifdef SWIFT_DEBUG_CHECKS
  c->cellID = pc->cellID;
#endif
363

364
  /* Copy the Multipole related information */
Matthieu Schaller's avatar
Matthieu Schaller committed
365
  if (with_gravity) {
366

367
    struct gravity_tensors *mp = c->grav.multipole;
368

369
370
371
372
373
374
375
376
377
    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;
378
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
379

380
381
382
383
  /* Number of new cells created. */
  int count = 1;

  /* Fill the progeny recursively, depth-first. */
384
  c->split = 0;
385
386
387
388
  for (int k = 0; k < 8; k++)
    if (pc->progeny[k] >= 0) {
      struct cell *temp;
      space_getcells(s, 1, &temp);
389
      temp->hydro.count = 0;
390
391
      temp->grav.count = 0;
      temp->stars.count = 0;
392
393
394
395
396
397
398
399
400
401
402
403
      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;
404
      temp->hydro.dx_max_part = 0.f;
405
      temp->hydro.dx_max_sort = 0.f;
Loic Hausammann's avatar
Loic Hausammann committed
406
      temp->stars.dx_max_part = 0.f;
Loic Hausammann's avatar
Loic Hausammann committed
407
      temp->stars.dx_max_sort = 0.f;
408
409
410
411
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
412
      count += cell_unpack(&pc[pc->progeny[k]], temp, s, with_gravity);
413
414
415
    }

  /* Return the total number of unpacked cells. */
416
  c->mpi.pcell_size = count;
417
418
419
420
421
422
423
424
  return count;

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

425
426
427
428
429
430
431
432
433
434
435
436
437
/**
 * @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. */
438
  c->mpi.tag = tags[0];
439
440
441
442
443
444
445
446
447
448
449

  /* 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
450
  if (c->mpi.pcell_size != count) error("Inconsistent tag and pcell count!");
451
452
453
454
455
456
457
458
459
460
461
#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
}

462
463
464
465
/**
 * @brief Pack the time information of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
466
 * @param pcells (output) The end-of-timestep information we pack into
467
468
469
 *
 * @return The number of packed cells.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
470
471
int cell_pack_end_step(struct cell *restrict c,
                       struct pcell_step *restrict pcells) {
472

473
474
#ifdef WITH_MPI

475
  /* Pack this cell's data. */
476
477
478
479
  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;
480
  pcells[0].stars.ti_end_min = c->stars.ti_end_min;
481
  pcells[0].hydro.dx_max_part = c->hydro.dx_max_part;
Loic Hausammann's avatar
Loic Hausammann committed
482
  pcells[0].stars.dx_max_part = c->stars.dx_max_part;
483

484
485
486
487
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
488
      count += cell_pack_end_step(c->progeny[k], &pcells[count]);
489
490
491
492
    }

  /* Return the number of packed values. */
  return count;
493
494
495
496
497

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

500
501
502
503
/**
 * @brief Unpack the time information of a given cell and its sub-cells.
 *
 * @param c The #cell
504
 * @param pcells The end-of-timestep information to unpack
505
506
507
 *
 * @return The number of cells created.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
508
509
int cell_unpack_end_step(struct cell *restrict c,
                         struct pcell_step *restrict pcells) {
510

511
512
#ifdef WITH_MPI

513
  /* Unpack this cell's data. */
514
515
516
517
  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;
518
  c->stars.ti_end_min = pcells[0].stars.ti_end_min;
519
  c->hydro.dx_max_part = pcells[0].hydro.dx_max_part;
Loic Hausammann's avatar
Loic Hausammann committed
520
  c->stars.dx_max_part = pcells[0].stars.dx_max_part;
521

522
523
524
525
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
526
      count += cell_unpack_end_step(c->progeny[k], &pcells[count]);
527
528
529
    }

  /* Return the number of packed values. */
530
  return count;
531
532
533
534
535

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

538
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
539
540
 * @brief Pack the multipole information of the given cell and all it's
 * sub-cells.
541
542
543
544
545
546
547
 *
 * @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
548
                         struct gravity_tensors *restrict pcells) {
549
550
551
552

#ifdef WITH_MPI

  /* Pack this cell's data. */
553
  pcells[0] = *c->grav.multipole;
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579

  /* 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
580
                           struct gravity_tensors *restrict pcells) {
581
582
583
584

#ifdef WITH_MPI

  /* Unpack this cell's data. */
585
  *c->grav.multipole = pcells[0];
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602

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

603
/**
604
 * @brief Lock a cell for access to its array of #part and hold its parents.
605
606
 *
 * @param c The #cell.
607
 * @return 0 on success, 1 on failure
608
 */
609
610
611
612
613
int cell_locktree(struct cell *c) {

  TIMER_TIC

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

  /* Did somebody hold this cell in the meantime? */
620
  if (c->hydro.hold) {
621
622

    /* Unlock this cell. */
623
    if (lock_unlock(&c->hydro.lock) != 0) error("Failed to unlock cell.");
624
625
626
627
628
629
630

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

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

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

    /* Increment the hold. */
638
    atomic_inc(&finger->hydro.hold);
639
640

    /* Unlock the cell. */
641
    if (lock_unlock(&finger->hydro.lock) != 0) error("Failed to unlock cell.");
642
643
644
645
646
647
648
649
650
651
652
653
  }

  /* 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
654
655
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
656
      atomic_dec(&finger2->hydro.hold);
657
658

    /* Unlock this cell. */
659
    if (lock_unlock(&c->hydro.lock) != 0) error("Failed to unlock cell.");
660
661
662
663
664
665
666

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

667
668
669
670
671
672
/**
 * @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
 */
673
674
675
676
677
int cell_glocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
678
  if (c->grav.phold || lock_trylock(&c->grav.plock) != 0) {
679
680
681
682
683
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
684
  if (c->grav.phold) {
685
686

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

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

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

    /* Lock this cell. */
699
    if (lock_trylock(&finger->grav.plock) != 0) break;
700
701

    /* Increment the hold. */
702
    atomic_inc(&finger->grav.phold);
703
704

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

  /* 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
718
719
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
720
      atomic_dec(&finger2->grav.phold);
721
722

    /* Unlock this cell. */
723
    if (lock_unlock(&c->grav.plock) != 0) error("Failed to unlock cell.");
724
725
726
727
728
729

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

731
732
733
734
735
736
737
738
739
740
741
/**
 * @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. */
742
  if (c->grav.mhold || lock_trylock(&c->grav.mlock) != 0) {
743
744
745
746
747
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
748
  if (c->grav.mhold) {
749
750

    /* Unlock this cell. */
751
    if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
752
753
754
755
756
757
758
759
760
761
762

    /* 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. */
763
    if (lock_trylock(&finger->grav.mlock) != 0) break;
764
765

    /* Increment the hold. */
766
    atomic_inc(&finger->grav.mhold);
767
768

    /* Unlock the cell. */
769
    if (lock_unlock(&finger->grav.mlock) != 0) error("Failed to unlock cell.");
770
771
772
773
774
775
776
777
778
779
780
781
782
783
  }

  /* 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)
784
      atomic_dec(&finger2->grav.mhold);
785
786

    /* Unlock this cell. */
787
    if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
788
789
790
791
792
793
794

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

795
796
797
798
799
800
801
802
803
804
805
/**
 * @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. */
806
  if (c->stars.hold || lock_trylock(&c->stars.lock) != 0) {
807
808
809
810
811
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
812
  if (c->stars.hold) {
813
814

    /* Unlock this cell. */
815
    if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
816
817
818
819
820
821
822
823
824
825
826

    /* 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. */
827
    if (lock_trylock(&finger->stars.lock) != 0) break;
828
829

    /* Increment the hold. */
830
    atomic_inc(&finger->stars.hold);
831
832

    /* Unlock the cell. */
833
    if (lock_unlock(&finger->stars.lock) != 0) error("Failed to unlock cell.");
834
835
836
837
838
839
840
841
842
843
844
845
846
847
  }

  /* 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)
848
      atomic_dec(&finger2->stars.hold);
849
850

    /* Unlock this cell. */
851
    if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
852
853
854
855
856
857
858

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

859
/**
860
 * @brief Unlock a cell's parents for access to #part array.
861
862
863
 *
 * @param c The #cell.
 */
864
865
866
867
868
void cell_unlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
872
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
873
    atomic_dec(&finger->hydro.hold);
874
875
876
877

  TIMER_TOC(timer_locktree);
}

878
879
880
881
882
/**
 * @brief Unlock a cell's parents for access to #gpart array.
 *
 * @param c The #cell.
 */
883
884
885
886
887
void cell_gunlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
891
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
892
    atomic_dec(&finger->grav.phold);
893
894
895
896

  TIMER_TOC(timer_locktree);
}

897
898
899
900
901
902
903
904
905
906
/**
 * @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. */
907
  if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
908
909
910

  /* Climb up the tree and unhold the parents. */
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
911
    atomic_dec(&finger->grav.mhold);
912
913
914
915

  TIMER_TOC(timer_locktree);
}

916
917
918
919
920
921
922
923
924
925
/**
 * @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. */
926
  if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
927
928
929

  /* Climb up the tree and unhold the parents. */
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
930
    atomic_dec(&finger->stars.hold);
931
932
933
934

  TIMER_TOC(timer_locktree);
}

935
936
937
938
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
939
 * @param parts_offset Offset of the cell parts array relative to the
940
 *        space's parts array, i.e. c->hydro.parts - s->parts.
941
 * @param sparts_offset Offset of the cell sparts array relative to the
942
943
 *        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)
944
 * entries, used for sorting indices.
945
946
947
 * @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)
948
 * entries, used for sorting indices for the gparts.
949
 */
950
951
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
                struct cell_buff *buff, struct cell_buff *sbuff,
952
                struct cell_buff *gbuff) {
953

954
955
  const int count = c->hydro.count, gcount = c->grav.count,
            scount = c->stars.count;
956
957
  struct part *parts = c->hydro.parts;
  struct xpart *xparts = c->hydro.xparts;
958
959
  struct gpart *gparts = c->grav.parts;
  struct spart *sparts = c->stars.parts;
960
961
962
963
964
965
  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];

966
967
968
#ifdef SWIFT_DEBUG_CHECKS
  /* Check that the buffs are OK. */
  for (int k = 0; k < count; k++) {
969
    if (buff[k].x[0] != parts[k].x[0] || buff[k].x[1] != parts[k].x[1] ||
970
        buff[k].x[2] != parts[k].x[2])
971
972
      error("Inconsistent buff contents.");
  }
973
974
975
976
977
978
979
980
981
982
  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.");
  }
983
#endif /* SWIFT_DEBUG_CHECKS */
984
985
986

  /* Fill the buffer with the indices. */
  for (int k = 0; k < count; k++) {
987
988
    const int bid = (buff[k].x[0] >= pivot[0]) * 4 +
                    (buff[k].x[1] >= pivot[1]) * 2 + (buff[k].x[2] >= pivot[2]);
989
    bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
990
    buff[k].ind = bid;
991
  }
992

993
994
995
996
997
  /* 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;
998
999
  }

1000
1001
1002
1003
  /* 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
1004
      int bid = buff[k].ind;
1005
1006
1007
      if (bid != bucket) {
        struct part part = parts[k];
        struct xpart xpart = xparts[k];
1008
        struct cell_buff temp_buff = buff[k];
1009
1010
        while (bid != bucket) {
          int j = bucket_offset[bid] + bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
1011
          while (buff[j].ind == bid) {
1012
1013
1014
            j++;
            bucket_count[bid]++;
          }
Pedro Gonnet's avatar
Pedro Gonnet committed
1015
1016
          memswap(&parts[j], &part, sizeof(struct part));
          memswap(&xparts[j], &xpart, sizeof(struct xpart));
1017
          memswap(&buff[j], &temp_buff, sizeof(struct cell_buff));
1018
1019
          if (parts[j].gpart)
            parts[j].gpart->id_or_neg_offset = -(j + parts_offset);
1020
          bid = temp_buff.ind;
1021
1022
1023
        }
        parts[k] = part;
        xparts[k] = xpart;
1024
        buff[k] = temp_buff;
1025
1026
        if (parts[k].gpart)
          parts[k].gpart->id_or_neg_offset = -(k + parts_offset);
1027
      }
1028
      bucket_count[bid]++;
1029
1030
1031
1032
    }
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
1033
  for (int k = 0; k < 8; k++) {
1034
    c->progeny[k]->hydro.count = bucket_count[k];
1035
    c->progeny[k]->hydro.count_total = c->progeny[k]->hydro.count;
1036
1037
    c->progeny[k]->hydro.parts = &c->hydro.parts[bucket_offset[k]];
    c->progeny[k]->hydro.xparts = &c->hydro.xparts[bucket_offset[k]];
1038
1039
  }

1040
#ifdef SWIFT_DEBUG_CHECKS
1041
1042
1043
1044
1045
1046
1047
1048
  /* 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);
  }

1049
  /* Verify that _all_ the parts have been assigned to a cell. */
1050
  for (int k = 1; k < 8; k++)
1051
1052
    if (&c->progeny[k - 1]->hydro.parts[c->progeny[k - 1]->hydro.count] !=
        c->progeny[k]->hydro.parts)
1053
      error("Particle sorting failed (internal consistency).");
1054
  if (c->progeny[0]->hydro.parts != c->hydro.parts)
1055
    error("Particle sorting failed (left edge).");
1056
1057
  if (&c->progeny[7]->hydro.parts[c->progeny[7]->hydro.count] !=
      &c->hydro.parts[count])
1058
    error("Particle sorting failed (right edge).");
1059
1060

  /* Verify a few sub-cells. */
1061
1062
1063
1064
  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])
1065
      error("Sorting failed (progeny=0).");
1066
1067
1068
1069
  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])
1070
      error("Sorting failed (progeny=1).");
1071
1072
1073
1074
  for (int k = 0; k < c->progeny[2]->hydro.count; k++)
    if (c->progeny[2]->hydro.parts[k].x[0] >= pivot[0] ||
        c->progeny[2]->hydro.parts[k].x[1] < pivot[1] ||
        c->progeny[2]->hydro.parts[k].x[2] >= pivot[2])
1075
      error("Sorting failed (progeny=2).");
1076
1077
1078
1079
  for (int k = 0; k < c->progeny[3]->hydro.count; k++)
    if (c->progeny[3]->hydro.parts[k].x[0] >= pivot[0] ||
        c->progeny[3]->hydro.parts[k].x[1] < pivot[1] ||
        c->progeny[3]->hydro.parts[k].x[2] < pivot[2])
1080
      error("Sorting failed (progeny=3).");
1081
1082
1083
1084
  for (int k = 0; k < c->progeny[4]->hydro.count; k++)
    if (c->progeny[4]->hydro.parts[k].x[0] < pivot[0] ||
        c->progeny[4]->hydro.parts[k].x[1] >= pivot[1] ||
        c->progeny[4]->hydro.parts[k].x[2] >= pivot[2])
1085
      error("Sorting failed (progeny=4).");
1086
1087
1088
1089
  for (int k = 0; k < c->progeny[5]->hydro.count; k++)
    if (c->progeny[5]->hydro.parts[k].x[0] < pivot[0] ||
        c->progeny[5]->hydro.parts[k].x[1] >= pivot[1] ||
        c->progeny[5]->hydro.parts[k].x[2] < pivot[2])
1090
      error("Sorting failed (progeny=5).");
1091
1092
1093
1094
  for (int k = 0; k < c->progeny[6]->hydro.count; k++)
    if (c->progeny[6]->hydro.parts[k].x[0] < pivot[0] ||
        c->progeny[6]->hydro.parts[k].x[1] < pivot[1] ||
        c->progeny[6]->hydro.parts[k].x[2] >= pivot[2])
1095
      error("Sorting failed (progeny=6).");
1096
1097
1098
1099
  for (int k = 0; k < c->progeny[7]->hydro.count; k++)
    if (c->progeny[7]->hydro.parts[k].x[0] < pivot[0] ||
        c->progeny[7]->hydro.parts[k].x[1] < pivot[1] ||
        c->progeny[7]->hydro.parts[k].x[2] < pivot[2])
1100
      error("Sorting failed (progeny=7).");
1101
#endif
1102