cell.c 156 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
int cell_link_parts(struct cell *c, struct part *parts) {

106
107
108
109
110
#ifdef SWIFT_DEBUG_CHECKS
  if(c->hydro.parts != NULL)
    error("Linking parts into a cell that was already linked");
#endif

111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
  c->hydro.parts = parts;

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

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


/**
 * @brief Link the cells recursively to the given #part array.
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */
int cell_link_foreign_parts(struct cell *c, struct part *parts) {
136
137
138
139
140

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

142
143
144
  /* Do we have a hydro task at this level? */
  if (c->hydro.density != NULL) {

145
146
147
148
149
150
151
    /* Recursively attach the parts */
    const int counts = cell_link_parts(c, parts);
#ifdef SWIFT_DEBUG_CHECKS
    if (counts != c->hydro.count)
      error("Something is wrong with the foreign counts");
#endif
    return counts;
152
  }
153

154
  /* Go deeper to find the level where the tasks are */
Pedro Gonnet's avatar
Pedro Gonnet committed
155
  if (c->split) {
156
    int count = 0;
Pedro Gonnet's avatar
Pedro Gonnet committed
157
    for (int k = 0; k < 8; k++) {
158
      if (c->progeny[k] != NULL) {
159
        count += cell_link_foreign_parts(c->progeny[k], &parts[count]);
160
      }
Pedro Gonnet's avatar
Pedro Gonnet committed
161
    }
162
    return count;
163
  } else {
164
    return 0;
165
  }
166
}
167

168
169
170
171
172
173
174
175
176
177
/**
 * @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) {

178
  c->grav.parts = gparts;
179
180
181
182
183
184
185
186
187
188
189

  /* 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. */
190
  return c->grav.count;
191
192
}

193
194
195
196
197
198
199
200
201
202
/**
 * @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) {

203
  c->stars.parts = sparts;
204
205
206
207
208
209
210
211
212
213
214

  /* 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. */
215
  return c->stars.count;
216
217
}

218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
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;
  }
}

247
248
249
250
251
252
/**
 * @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.
253
254
 * @param with_gravity Are we running with gravity and hence need
 *      to exchange multipoles?
255
256
257
 *
 * @return The number of packed cells.
 */
258
int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
Matthieu Schaller's avatar
Matthieu Schaller committed
259
              const int with_gravity) {
260

261
262
#ifdef WITH_MPI

263
  /* Start by packing the data of the current cell. */
264
265
266
267
268
  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;
269
  pc->stars.ti_end_min = c->stars.ti_end_min;
270
271
  pc->hydro.ti_old_part = c->hydro.ti_old_part;
  pc->grav.ti_old_part = c->grav.ti_old_part;
272
273
  pc->grav.ti_old_multipole = c->grav.ti_old_multipole;
  pc->hydro.count = c->hydro.count;
274
275
  pc->grav.count = c->grav.count;
  pc->stars.count = c->stars.count;
276
  pc->maxdepth = c->maxdepth;
277

278
  /* Copy the Multipole related information */
Matthieu Schaller's avatar
Matthieu Schaller committed
279
  if (with_gravity) {
280
    const struct gravity_tensors *mp = c->grav.multipole;
281

282
283
284
285
286
287
288
289
290
    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;
291
292
  }

293
294
295
#ifdef SWIFT_DEBUG_CHECKS
  pc->cellID = c->cellID;
#endif
296
297

  /* Fill in the progeny, depth-first recursion. */
Pedro Gonnet's avatar
Pedro Gonnet committed
298
299
  int count = 1;
  for (int k = 0; k < 8; k++)
300
301
    if (c->progeny[k] != NULL) {
      pc->progeny[k] = count;
302
      count += cell_pack(c->progeny[k], &pc[count], with_gravity);
303
    } else {
304
      pc->progeny[k] = -1;
305
    }
306
307

  /* Return the number of packed cells used. */
308
  c->mpi.pcell_size = count;
309
  return count;
310
311
312
313
314

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

317
318
319
320
321
322
323
324
325
326
327
328
329
/**
 * @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. */
330
  tags[0] = c->mpi.tag;
331
332
333
334
335
336
337
338

  /* 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
339
  if (c->mpi.pcell_size != count) error("Inconsistent tag and pcell count!");
340
341
342
343
344
345
346
347
348
349
350
#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
}

351
352
353
354
355
356
/**
 * @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.
357
358
 * @param with_gravity Are we running with gravity and hence need
 *      to exchange multipoles?
359
360
361
 *
 * @return The number of cells created.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
362
int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
363
                struct space *restrict s, const int with_gravity) {
364
365
366
367

#ifdef WITH_MPI

  /* Unpack the current pcell. */
368
369
370
371
372
  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;
373
  c->stars.ti_end_min = pc->stars.ti_end_min;
374
375
  c->hydro.ti_old_part = pc->hydro.ti_old_part;
  c->grav.ti_old_part = pc->grav.ti_old_part;
376
377
  c->grav.ti_old_multipole = pc->grav.ti_old_multipole;
  c->hydro.count = pc->hydro.count;
378
379
  c->grav.count = pc->grav.count;
  c->stars.count = pc->stars.count;
380
381
  c->maxdepth = pc->maxdepth;

382
383
384
#ifdef SWIFT_DEBUG_CHECKS
  c->cellID = pc->cellID;
#endif
385

386
  /* Copy the Multipole related information */
Matthieu Schaller's avatar
Matthieu Schaller committed
387
  if (with_gravity) {
388

389
    struct gravity_tensors *mp = c->grav.multipole;
390

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

402
403
404
405
  /* Number of new cells created. */
  int count = 1;

  /* Fill the progeny recursively, depth-first. */
406
  c->split = 0;
407
408
409
410
  for (int k = 0; k < 8; k++)
    if (pc->progeny[k] >= 0) {
      struct cell *temp;
      space_getcells(s, 1, &temp);
411
      temp->hydro.count = 0;
412
413
      temp->grav.count = 0;
      temp->stars.count = 0;
414
415
416
417
418
419
420
421
422
423
424
425
      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;
426
      temp->hydro.dx_max_part = 0.f;
427
      temp->hydro.dx_max_sort = 0.f;
Loic Hausammann's avatar
Loic Hausammann committed
428
      temp->stars.dx_max_part = 0.f;
Loic Hausammann's avatar
Loic Hausammann committed
429
      temp->stars.dx_max_sort = 0.f;
430
431
432
433
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
434
      count += cell_unpack(&pc[pc->progeny[k]], temp, s, with_gravity);
435
436
437
    }

  /* Return the total number of unpacked cells. */
438
  c->mpi.pcell_size = count;
439
440
441
442
443
444
445
446
  return count;

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

447
448
449
450
451
452
453
454
455
456
457
458
459
/**
 * @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. */
460
  c->mpi.tag = tags[0];
461
462
463
464
465
466
467
468
469
470
471

  /* 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
472
  if (c->mpi.pcell_size != count) error("Inconsistent tag and pcell count!");
473
474
475
476
477
478
479
480
481
482
483
#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
}

484
485
486
487
/**
 * @brief Pack the time information of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
488
 * @param pcells (output) The end-of-timestep information we pack into
489
490
491
 *
 * @return The number of packed cells.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
492
493
int cell_pack_end_step(struct cell *restrict c,
                       struct pcell_step *restrict pcells) {
494

495
496
#ifdef WITH_MPI

497
  /* Pack this cell's data. */
498
499
500
501
  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;
502
  pcells[0].stars.ti_end_min = c->stars.ti_end_min;
503
  pcells[0].hydro.dx_max_part = c->hydro.dx_max_part;
Loic Hausammann's avatar
Loic Hausammann committed
504
  pcells[0].stars.dx_max_part = c->stars.dx_max_part;
505

506
507
508
509
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
510
      count += cell_pack_end_step(c->progeny[k], &pcells[count]);
511
512
513
514
    }

  /* Return the number of packed values. */
  return count;
515
516
517
518
519

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

522
523
524
525
/**
 * @brief Unpack the time information of a given cell and its sub-cells.
 *
 * @param c The #cell
526
 * @param pcells The end-of-timestep information to unpack
527
528
529
 *
 * @return The number of cells created.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
530
531
int cell_unpack_end_step(struct cell *restrict c,
                         struct pcell_step *restrict pcells) {
532

533
534
#ifdef WITH_MPI

535
  /* Unpack this cell's data. */
536
537
538
539
  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;
540
  c->stars.ti_end_min = pcells[0].stars.ti_end_min;
541
  c->hydro.dx_max_part = pcells[0].hydro.dx_max_part;
Loic Hausammann's avatar
Loic Hausammann committed
542
  c->stars.dx_max_part = pcells[0].stars.dx_max_part;
543

544
545
546
547
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
548
      count += cell_unpack_end_step(c->progeny[k], &pcells[count]);
549
550
551
    }

  /* Return the number of packed values. */
552
  return count;
553
554
555
556
557

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

560
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
561
562
 * @brief Pack the multipole information of the given cell and all it's
 * sub-cells.
563
564
565
566
567
568
569
 *
 * @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
570
                         struct gravity_tensors *restrict pcells) {
571
572
573
574

#ifdef WITH_MPI

  /* Pack this cell's data. */
575
  pcells[0] = *c->grav.multipole;
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601

  /* 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
602
                           struct gravity_tensors *restrict pcells) {
603
604
605
606

#ifdef WITH_MPI

  /* Unpack this cell's data. */
607
  *c->grav.multipole = pcells[0];
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624

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

625
/**
626
 * @brief Lock a cell for access to its array of #part and hold its parents.
627
628
 *
 * @param c The #cell.
629
 * @return 0 on success, 1 on failure
630
 */
631
632
633
634
635
int cell_locktree(struct cell *c) {

  TIMER_TIC

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

  /* Did somebody hold this cell in the meantime? */
642
  if (c->hydro.hold) {
643
644

    /* Unlock this cell. */
645
    if (lock_unlock(&c->hydro.lock) != 0) error("Failed to unlock cell.");
646
647
648
649
650
651
652

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

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

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

    /* Increment the hold. */
660
    atomic_inc(&finger->hydro.hold);
661
662

    /* Unlock the cell. */
663
    if (lock_unlock(&finger->hydro.lock) != 0) error("Failed to unlock cell.");
664
665
666
667
668
669
670
671
672
673
674
675
  }

  /* 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
676
677
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
678
      atomic_dec(&finger2->hydro.hold);
679
680

    /* Unlock this cell. */
681
    if (lock_unlock(&c->hydro.lock) != 0) error("Failed to unlock cell.");
682
683
684
685
686
687
688

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

689
690
691
692
693
694
/**
 * @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
 */
695
696
697
698
699
int cell_glocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
700
  if (c->grav.phold || lock_trylock(&c->grav.plock) != 0) {
701
702
703
704
705
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
706
  if (c->grav.phold) {
707
708

    /* Unlock this cell. */
709
    if (lock_unlock(&c->grav.plock) != 0) error("Failed to unlock cell.");
710
711
712
713
714
715
716

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

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

    /* Lock this cell. */
721
    if (lock_trylock(&finger->grav.plock) != 0) break;
722
723

    /* Increment the hold. */
724
    atomic_inc(&finger->grav.phold);
725
726

    /* Unlock the cell. */
727
    if (lock_unlock(&finger->grav.plock) != 0) error("Failed to unlock cell.");
728
729
730
731
732
733
734
735
736
737
738
739
  }

  /* 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
740
741
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
742
      atomic_dec(&finger2->grav.phold);
743
744

    /* Unlock this cell. */
745
    if (lock_unlock(&c->grav.plock) != 0) error("Failed to unlock cell.");
746
747
748
749
750
751

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

753
754
755
756
757
758
759
760
761
762
763
/**
 * @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. */
764
  if (c->grav.mhold || lock_trylock(&c->grav.mlock) != 0) {
765
766
767
768
769
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
770
  if (c->grav.mhold) {
771
772

    /* Unlock this cell. */
773
    if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
774
775
776
777
778
779
780
781
782
783
784

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

    /* Increment the hold. */
788
    atomic_inc(&finger->grav.mhold);
789
790

    /* Unlock the cell. */
791
    if (lock_unlock(&finger->grav.mlock) != 0) error("Failed to unlock cell.");
792
793
794
795
796
797
798
799
800
801
802
803
804
805
  }

  /* 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)
806
      atomic_dec(&finger2->grav.mhold);
807
808

    /* Unlock this cell. */
809
    if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
810
811
812
813
814
815
816

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

817
818
819
820
821
822
823
824
825
826
827
/**
 * @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. */
828
  if (c->stars.hold || lock_trylock(&c->stars.lock) != 0) {
829
830
831
832
833
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
834
  if (c->stars.hold) {
835
836

    /* Unlock this cell. */
837
    if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
838
839
840
841
842
843
844
845
846
847
848

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

    /* Increment the hold. */
852
    atomic_inc(&finger->stars.hold);
853
854

    /* Unlock the cell. */
855
    if (lock_unlock(&finger->stars.lock) != 0) error("Failed to unlock cell.");
856
857
858
859
860
861
862
863
864
865
866
867
868
869
  }

  /* 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)
870
      atomic_dec(&finger2->stars.hold);
871
872

    /* Unlock this cell. */
873
    if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
874
875
876
877
878
879
880

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

881
/**
882
 * @brief Unlock a cell's parents for access to #part array.
883
884
885
 *
 * @param c The #cell.
 */
886
887
888
889
890
void cell_unlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
894
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
895
    atomic_dec(&finger->hydro.hold);
896
897
898
899

  TIMER_TOC(timer_locktree);
}

900
901
902
903
904
/**
 * @brief Unlock a cell's parents for access to #gpart array.
 *
 * @param c The #cell.
 */
905
906
907
908
909
void cell_gunlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
913
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
914
    atomic_dec(&finger->grav.phold);
915
916
917
918

  TIMER_TOC(timer_locktree);
}

919
920
921
922
923
924
925
926
927
928
/**
 * @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. */
929
  if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
930
931
932

  /* Climb up the tree and unhold the parents. */
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
933
    atomic_dec(&finger->grav.mhold);
934
935
936
937

  TIMER_TOC(timer_locktree);
}

938
939
940
941
942
943
944
945
946
947
/**
 * @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. */
948
  if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
949
950
951

  /* Climb up the tree and unhold the parents. */
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
952
    atomic_dec(&finger->stars.hold);
953
954
955
956

  TIMER_TOC(timer_locktree);
}

957
958
959
960
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
961
 * @param parts_offset Offset of the cell parts array relative to the
962
 *        space's parts array, i.e. c->hydro.parts - s->parts.
963
 * @param sparts_offset Offset of the cell sparts array relative to the
964
965
 *        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)
966
 * entries, used for sorting indices.
967
968
969
 * @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)
970
 * entries, used for sorting indices for the gparts.
971
 */
972
973
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
                struct cell_buff *buff, struct cell_buff *sbuff,
974
                struct cell_buff *gbuff) {
975

976
977
  const int count = c->hydro.count, gcount = c->grav.count,
            scount = c->stars.count;
978
979
  struct part *parts = c->hydro.parts;
  struct xpart *xparts = c->hydro.xparts;
980
981
  struct gpart *gparts = c->grav.parts;
  struct spart *sparts = c->stars.parts;
982
983
984
985
986
987
  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];

988
989
990
#ifdef SWIFT_DEBUG_CHECKS
  /* Check that the buffs are OK. */
  for (int k = 0; k < count; k++) {
991
    if (buff[k].x[0] != parts[k].x[0] || buff[k].x[1] != parts[k].x[1] ||
992
        buff[k].x[2] != parts[k].x[2])
993
994
      error("Inconsistent buff contents.");
  }
995
996
997
998
999
1000
1001
1002
1003
1004
  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.");
  }
1005
#endif /* SWIFT_DEBUG_CHECKS */
1006
1007
1008

  /* Fill the buffer with the indices. */
  for (int k = 0; k < count; k++) {
1009
1010
    const int bid = (buff[k].x[0] >= pivot[0]) * 4 +
                    (buff[k].x[1] >= pivot[1]) * 2 + (buff[k].x[2] >= pivot[2]);
1011
    bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
1012
    buff[k].ind = bid;
1013
  }
1014

1015
1016
1017
1018
1019
  /* 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;
1020
1021
  }

1022
1023
1024
1025
  /* 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
1026
      int bid = buff[k].ind;
1027
1028
1029
      if (bid != bucket) {
        struct part part = parts[k];
        struct xpart xpart = xparts[k];
1030
        struct cell_buff temp_buff = buff[k];
1031
1032
        while (bid != bucket) {
          int j = bucket_offset[bid] + bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
1033
          while (buff[j].ind == bid) {
1034
1035
1036
            j++;
            bucket_count[bid]++;
          }
Pedro Gonnet's avatar
Pedro Gonnet committed
1037
1038
          memswap(&parts[j], &part, sizeof(struct part));
          memswap(&xparts[j], &xpart, sizeof(struct xpart));
1039
          memswap(&buff[j], &temp_buff, sizeof(struct cell_buff));
1040
1041
          if (parts[j].gpart)
            parts[j].gpart->id_or_neg_offset = -(j + parts_offset);
1042
          bid = temp_buff.ind;
1043
1044
1045
        }
        parts[k] = part;
        xparts[k] = xpart;
1046
        buff[k] = temp_buff;
1047
1048
        if (parts[k].gpart)
          parts[k].gpart->id_or_neg_offset = -(k + parts_offset);
1049
      }
1050
      bucket_count[bid]++;
1051
1052
1053
1054
    }
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
1055
  for (int k = 0; k < 8; k++) {
1056
    c->progeny[k]->hydro.count = bucket_count[k];
1057
    c->progeny[k]->hydro.count_total = c->progeny[k]->hydro.count;
1058
1059
    c->progeny[k]->hydro.parts = &c->hydro.parts[bucket_offset[k]];
    c->progeny[k]->hydro.xparts = &c->hydro.xparts[bucket_offset[k]];
1060
1061
  }

1062
#ifdef SWIFT_DEBUG_CHECKS
1063
1064
1065
1066
1067
1068
1069
1070
  /* 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);
  }

1071
  /* Verify that _all_ the parts have been assigned to a cell. */
1072
  for (int k = 1; k < 8; k++)
1073
1074
    if (&c->progeny[k - 1]->hydro.parts[c->progeny[k - 1]->hydro.count] !=
        c->progeny[k]->hydro.parts)
1075
      error("Particle sorting failed (internal consistency).");
1076
  if (c->progeny[0]->hydro.parts != c->hydro.parts)
1077
    error("Particle sorting failed (left edge).");
1078
1079
  if (&c->progeny[7]->hydro.parts[c->progeny[7]->hydro.count] !=
      &c->hydro.parts[count])
1080
    error("Particle sorting failed (right edge).");
1081
1082

  /* Verify a few sub-cells. */
1083
1084
1085
1086
  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])
1087
      error("Sorting failed (progeny=0).");
1088
1089
1090
1091
  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])
1092
      error("Sorting failed (progeny=1).");
1093
1094
1095
1096
  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])
1097
      error("Sorting failed (progeny=2).");
1098
1099
1100
1101
  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])
1102
      error("Sorting failed (progeny=3).");
1103
1104
1105
1106
  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])
1107
      error("Sorting failed (progeny=4).");
1108
1109
1110
1111
  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])
1112