cell.c 133 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])
1001
      error("Sorting failed (progeny=1).");
1002
1003
1004
1005
  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])
1006
      error("Sorting failed (progeny=2).");
1007
1008
1009
1010
  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])
1011
      error("Sorting failed (progeny=3).");
1012
1013
1014
1015
  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])
1016
      error("Sorting failed (progeny=4).");
1017
1018
1019
1020
  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])
1021
      error("Sorting failed (progeny=5).");
1022
1023
1024
1025
  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])
1026
      error("Sorting failed (progeny=6).");
1027
1028
1029
1030
  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])
1031
      error("Sorting failed (progeny=7).");
1032
#endif
1033

1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
  /* Now do the same song and dance for the sparts. */
  for (int k = 0; k < 8; k++) bucket_count[k] = 0;

  /* Fill the buffer with the indices. */
  for (int k = 0; k < scount; k++) {
    const int bid = (sbuff[k].x[0] > pivot[0]) * 4 +
                    (sbuff[k].x[1] > pivot[1]) * 2 + (sbuff[k].x[2] > pivot[2]);
    bucket_count[bid]++;
    sbuff[k].ind = bid;
  }

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

  /* 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++) {
      int bid = sbuff[k].ind;
      if (bid != bucket) {
        struct spart spart = sparts[k];
        struct cell_buff temp_buff = sbuff[k];
        while (bid != bucket) {
          int j = bucket_offset[bid] + bucket_count[bid]++;
          while (sbuff[j].ind == bid) {
            j++;
            bucket_count[bid]++;
          }
          memswap(&sparts[j], &spart, sizeof(struct spart));
          memswap(&sbuff[j], &temp_buff, sizeof(struct cell_buff));
1068
1069
          if (sparts[j].gpart)
            sparts[j].gpart->id_or_neg_offset = -(j + sparts_offset);
1070
1071
1072
1073
          bid = temp_buff.ind;
        }
        sparts[k] = spart;
        sbuff[k] = temp_buff;
1074
1075
        if (sparts[k].gpart)
          sparts[k].gpart->id_or_neg_offset = -(k + sparts_offset);
1076
1077
1078
1079
1080
1081
1082
      }
      bucket_count[bid]++;
    }
  }

  /* Store the counts and offsets. */
  for (int k = 0; k < 8; k++) {
1083
1084
    c->progeny[k]->stars.count = bucket_count[k];
    c->progeny[k]->stars.parts = &c->stars.parts[bucket_offset[k]];
1085
1086
1087
  }

  /* Finally, do the same song and dance for the gparts. */
1088
1089
1090
1091
  for (int k = 0; k < 8; k++) bucket_count[k] = 0;

  /* Fill the buffer with the indices. */
  for (int k = 0; k < gcount; k++) {
1092
1093
    const int bid = (gbuff[k].x[0] > pivot[0]) * 4 +
                    (gbuff[k].x[1] > pivot[1]) * 2 + (gbuff[k].x[2] > pivot[2]);
1094
    bucket_count[bid]++;
1095
    gbuff[k].ind = bid;
1096
  }
1097
1098
1099
1100
1101
1102

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