cell.c 168 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
/**
 * @brief Get the size of the cell subtree.
 *
 * @param c The #cell.
 */
77
int cell_getsize(struct cell *c) {
78

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

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

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

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

101
#ifdef SWIFT_DEBUG_CHECKS
102
103
104
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");

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

109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
  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;
}

/**
125
 * @brief Link the cells recursively to the given #gpart array.
126
127
 *
 * @param c The #cell.
128
 * @param gparts The #gpart array.
129
130
131
 *
 * @return The number of particles linked.
 */
132
int cell_link_gparts(struct cell *c, struct gpart *gparts) {
133
134
135
136
137

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

138
  if (c->grav.parts != NULL)
139
    error("Linking gparts into a cell that was already linked");
140
#endif
141

142
  c->grav.parts = gparts;
143
144
145
146
147
148
149
150
151
152
153

  /* 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. */
154
  return c->grav.count;
155
156
}

157
158
159
160
161
162
163
164
165
166
/**
 * @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) {

167
168
169
170
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");

171
  if (c->stars.parts != NULL)
172
173
174
    error("Linking sparts into a cell that was already linked");
#endif

175
  c->stars.parts = sparts;
176
177
178
179
180
181
182
183
184
185
186

  /* 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. */
187
  return c->stars.count;
188
189
}

190
191
192
193
194
195
196
197
198
199
/**
 * @brief Recurse down foreign cells until reaching one with hydro
 * tasks; then trigger the linking of the #part array from that
 * level.
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */
200
201
int cell_link_foreign_parts(struct cell *c, struct part *parts) {

202
203
#ifdef WITH_MPI

204
205
206
207
208
209
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");
#endif

  /* Do we have a hydro task at this level? */
210
  if (c->mpi.hydro.recv_xv != NULL) {
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228

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

  /* Go deeper to find the level where the tasks are */
  if (c->split) {
    int count = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) {
        count += cell_link_foreign_parts(c->progeny[k], &parts[count]);
      }
    }
229
230
231
    return count;
  } else {
    return 0;
232
  }
233
234
235
236

#else
  error("Calling linking of foregin particles in non-MPI mode.");
#endif
237
238
}

239
240
241
242
243
244
245
246
247
248
/**
 * @brief Recurse down foreign cells until reaching one with gravity
 * tasks; then trigger the linking of the #gpart array from that
 * level.
 *
 * @param c The #cell.
 * @param gparts The #gpart array.
 *
 * @return The number of particles linked.
 */
249
250
int cell_link_foreign_gparts(struct cell *c, struct gpart *gparts) {

251
252
#ifdef WITH_MPI

253
254
255
256
257
258
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");
#endif

  /* Do we have a hydro task at this level? */
259
  if (c->mpi.grav.recv != NULL) {
260

261
    /* Recursively attach the gparts */
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
    const int counts = cell_link_gparts(c, gparts);
#ifdef SWIFT_DEBUG_CHECKS
    if (counts != c->grav.count)
      error("Something is wrong with the foreign counts");
#endif
    return counts;
  }

  /* Go deeper to find the level where the tasks are */
  if (c->split) {
    int count = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) {
        count += cell_link_foreign_gparts(c->progeny[k], &gparts[count]);
      }
    }
278
279
280
    return count;
  } else {
    return 0;
281
  }
282
283
284
285

#else
  error("Calling linking of foregin particles in non-MPI mode.");
#endif
286
287
}

288
289
290
291
292
293
294
295
/**
 * @brief Recursively count the number of #part in foreign cells that
 * are in cells with hydro-related tasks.
 *
 * @param c The #cell.
 *
 * @return The number of particles linked.
 */
296
297
int cell_count_parts_for_tasks(const struct cell *c) {

298
299
#ifdef WITH_MPI

300
301
302
303
304
305
#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? */
306
  if (c->mpi.hydro.recv_xv != NULL) {
307
308
309
310
311
312
313
314
315
316
    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]);
      }
    }
317
318
319
    return count;
  } else {
    return 0;
320
  }
321
322
323
324

#else
  error("Calling linking of foregin particles in non-MPI mode.");
#endif
325
326
}

327
328
329
330
331
332
333
334
/**
 * @brief Recursively count the number of #gpart in foreign cells that
 * are in cells with gravity-related tasks.
 *
 * @param c The #cell.
 *
 * @return The number of particles linked.
 */
335
336
int cell_count_gparts_for_tasks(const struct cell *c) {

337
338
#ifdef WITH_MPI

339
340
341
342
343
344
#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? */
345
  if (c->mpi.grav.recv != NULL) {
346
347
348
349
350
351
352
353
354
355
    return c->grav.count;
  }

  if (c->split) {
    int count = 0;
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
        count += cell_count_gparts_for_tasks(c->progeny[k]);
      }
    }
356
357
358
    return count;
  } else {
    return 0;
359
  }
360
361
362
363

#else
  error("Calling linking of foregin particles in non-MPI mode.");
#endif
364
365
}

366
367
368
369
370
371
/**
 * @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.
372
373
 * @param with_gravity Are we running with gravity and hence need
 *      to exchange multipoles?
374
375
376
 *
 * @return The number of packed cells.
 */
377
int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
Matthieu Schaller's avatar
Matthieu Schaller committed
378
              const int with_gravity) {
379

380
381
#ifdef WITH_MPI

382
  /* Start by packing the data of the current cell. */
383
  pc->hydro.h_max = c->hydro.h_max;
384
  pc->stars.h_max = c->stars.h_max;
385
386
387
388
  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;
389
  pc->stars.ti_end_min = c->stars.ti_end_min;
390
  pc->stars.ti_end_max = c->stars.ti_end_max;
391
392
  pc->hydro.ti_old_part = c->hydro.ti_old_part;
  pc->grav.ti_old_part = c->grav.ti_old_part;
393
  pc->grav.ti_old_multipole = c->grav.ti_old_multipole;
394
  pc->stars.ti_old_part = c->stars.ti_old_part;
395
  pc->hydro.count = c->hydro.count;
396
397
  pc->grav.count = c->grav.count;
  pc->stars.count = c->stars.count;
398
  pc->maxdepth = c->maxdepth;
399

400
  /* Copy the Multipole related information */
Matthieu Schaller's avatar
Matthieu Schaller committed
401
  if (with_gravity) {
402
    const struct gravity_tensors *mp = c->grav.multipole;
403

404
405
406
407
408
409
410
411
412
    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;
413
414
  }

415
416
417
#ifdef SWIFT_DEBUG_CHECKS
  pc->cellID = c->cellID;
#endif
418
419

  /* Fill in the progeny, depth-first recursion. */
Pedro Gonnet's avatar
Pedro Gonnet committed
420
421
  int count = 1;
  for (int k = 0; k < 8; k++)
422
423
    if (c->progeny[k] != NULL) {
      pc->progeny[k] = count;
424
      count += cell_pack(c->progeny[k], &pc[count], with_gravity);
425
    } else {
426
      pc->progeny[k] = -1;
427
    }
428
429

  /* Return the number of packed cells used. */
430
  c->mpi.pcell_size = count;
431
  return count;
432
433
434
435
436

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

439
440
441
442
443
444
445
446
447
448
449
450
451
/**
 * @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. */
452
  tags[0] = c->mpi.tag;
453
454
455
456
457
458
459
460

  /* 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
461
  if (c->mpi.pcell_size != count) error("Inconsistent tag and pcell count!");
462
463
464
465
466
467
468
469
470
471
472
#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
}

473
474
475
476
477
478
/**
 * @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.
479
480
 * @param with_gravity Are we running with gravity and hence need
 *      to exchange multipoles?
481
482
483
 *
 * @return The number of cells created.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
484
int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
485
                struct space *restrict s, const int with_gravity) {
486
487
488
489

#ifdef WITH_MPI

  /* Unpack the current pcell. */
490
  c->hydro.h_max = pc->hydro.h_max;
491
  c->stars.h_max = pc->stars.h_max;
492
493
494
495
  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;
496
  c->stars.ti_end_min = pc->stars.ti_end_min;
497
  c->stars.ti_end_max = pc->stars.ti_end_max;
498
499
  c->hydro.ti_old_part = pc->hydro.ti_old_part;
  c->grav.ti_old_part = pc->grav.ti_old_part;
500
  c->grav.ti_old_multipole = pc->grav.ti_old_multipole;
501
  c->stars.ti_old_part = pc->stars.ti_old_part;
502
  c->hydro.count = pc->hydro.count;
503
504
  c->grav.count = pc->grav.count;
  c->stars.count = pc->stars.count;
505
506
  c->maxdepth = pc->maxdepth;

507
508
509
#ifdef SWIFT_DEBUG_CHECKS
  c->cellID = pc->cellID;
#endif
510

511
  /* Copy the Multipole related information */
Matthieu Schaller's avatar
Matthieu Schaller committed
512
  if (with_gravity) {
513

514
    struct gravity_tensors *mp = c->grav.multipole;
515

516
517
518
519
520
521
522
523
524
    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;
525
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
526

527
528
529
530
  /* Number of new cells created. */
  int count = 1;

  /* Fill the progeny recursively, depth-first. */
531
  c->split = 0;
532
533
534
535
  for (int k = 0; k < 8; k++)
    if (pc->progeny[k] >= 0) {
      struct cell *temp;
      space_getcells(s, 1, &temp);
536
      temp->hydro.count = 0;
537
538
      temp->grav.count = 0;
      temp->stars.count = 0;
539
540
541
542
543
544
545
546
547
548
549
550
      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;
551
      temp->hydro.dx_max_part = 0.f;
552
      temp->hydro.dx_max_sort = 0.f;
Loic Hausammann's avatar
Loic Hausammann committed
553
      temp->stars.dx_max_part = 0.f;
Loic Hausammann's avatar
Loic Hausammann committed
554
      temp->stars.dx_max_sort = 0.f;
555
556
557
558
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
559
      count += cell_unpack(&pc[pc->progeny[k]], temp, s, with_gravity);
560
561
562
    }

  /* Return the total number of unpacked cells. */
563
  c->mpi.pcell_size = count;
564
565
566
567
568
569
570
571
  return count;

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

572
573
574
575
576
577
578
579
580
581
582
583
584
/**
 * @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. */
585
  c->mpi.tag = tags[0];
586
587
588
589
590
591
592
593
594
595
596

  /* 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
597
  if (c->mpi.pcell_size != count) error("Inconsistent tag and pcell count!");
598
599
600
601
602
603
604
605
606
607
608
#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
}

609
610
611
612
/**
 * @brief Pack the time information of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
613
 * @param pcells (output) The end-of-timestep information we pack into
614
615
616
 *
 * @return The number of packed cells.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
617
618
int cell_pack_end_step(struct cell *restrict c,
                       struct pcell_step *restrict pcells) {
619

620
621
#ifdef WITH_MPI

622
  /* Pack this cell's data. */
623
624
625
626
  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;
627
  pcells[0].stars.ti_end_min = c->stars.ti_end_min;
628
  pcells[0].stars.ti_end_max = c->stars.ti_end_max;
629
  pcells[0].hydro.dx_max_part = c->hydro.dx_max_part;
Loic Hausammann's avatar
Loic Hausammann committed
630
  pcells[0].stars.dx_max_part = c->stars.dx_max_part;
631

632
633
634
635
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
636
      count += cell_pack_end_step(c->progeny[k], &pcells[count]);
637
638
639
640
    }

  /* Return the number of packed values. */
  return count;
641
642
643
644
645

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

648
649
650
651
/**
 * @brief Unpack the time information of a given cell and its sub-cells.
 *
 * @param c The #cell
652
 * @param pcells The end-of-timestep information to unpack
653
654
655
 *
 * @return The number of cells created.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
656
657
int cell_unpack_end_step(struct cell *restrict c,
                         struct pcell_step *restrict pcells) {
658

659
660
#ifdef WITH_MPI

661
  /* Unpack this cell's data. */
662
663
664
665
  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;
666
  c->stars.ti_end_min = pcells[0].stars.ti_end_min;
667
  c->stars.ti_end_max = pcells[0].stars.ti_end_max;
668
  c->hydro.dx_max_part = pcells[0].hydro.dx_max_part;
Loic Hausammann's avatar
Loic Hausammann committed
669
  c->stars.dx_max_part = pcells[0].stars.dx_max_part;
670

671
672
673
674
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
675
      count += cell_unpack_end_step(c->progeny[k], &pcells[count]);
676
677
678
    }

  /* Return the number of packed values. */
679
  return count;
680
681
682
683
684

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

687
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
688
689
 * @brief Pack the multipole information of the given cell and all it's
 * sub-cells.
690
691
692
693
694
695
696
 *
 * @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
697
                         struct gravity_tensors *restrict pcells) {
698
699
700
701

#ifdef WITH_MPI

  /* Pack this cell's data. */
702
  pcells[0] = *c->grav.multipole;
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728

  /* 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
729
                           struct gravity_tensors *restrict pcells) {
730
731
732
733

#ifdef WITH_MPI

  /* Unpack this cell's data. */
734
  *c->grav.multipole = pcells[0];
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751

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

752
/**
753
 * @brief Lock a cell for access to its array of #part and hold its parents.
754
755
 *
 * @param c The #cell.
756
 * @return 0 on success, 1 on failure
757
 */
758
759
760
761
762
int cell_locktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
763
  if (c->hydro.hold || lock_trylock(&c->hydro.lock) != 0) {
764
765
766
767
768
    TIMER_TOC(timer_locktree);
    return 1;
  }

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

    /* Unlock this cell. */
772
    if (lock_unlock(&c->hydro.lock) != 0) error("Failed to unlock cell.");
773
774
775
776
777
778
779

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

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

    /* Lock this cell. */
784
    if (lock_trylock(&finger->hydro.lock) != 0) break;
785
786

    /* Increment the hold. */
787
    atomic_inc(&finger->hydro.hold);
788
789

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

  /* 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
803
804
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
805
      atomic_dec(&finger2->hydro.hold);
806
807

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

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

816
817
818
819
820
821
/**
 * @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
 */
822
823
824
825
826
int cell_glocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
827
  if (c->grav.phold || lock_trylock(&c->grav.plock) != 0) {
828
829
830
831
832
    TIMER_TOC(timer_locktree);
    return 1;
  }

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

    /* Unlock this cell. */
836
    if (lock_unlock(&c->grav.plock) != 0) error("Failed to unlock cell.");
837
838
839
840
841
842
843

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

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

    /* Lock this cell. */
848
    if (lock_trylock(&finger->grav.plock) != 0) break;
849
850

    /* Increment the hold. */
851
    atomic_inc(&finger->grav.phold);
852
853

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

  /* 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
867
868
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
869
      atomic_dec(&finger2->grav.phold);
870
871

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

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

880
881
882
883
884
885
886
887
888
889
890
/**
 * @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. */
891
  if (c->grav.mhold || lock_trylock(&c->grav.mlock) != 0) {
892
893
894
895
896
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
897
  if (c->grav.mhold) {
898
899

    /* Unlock this cell. */
900
    if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
901
902
903
904
905
906
907
908
909
910
911

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

    /* Increment the hold. */
915
    atomic_inc(&finger->grav.mhold);
916
917

    /* Unlock the cell. */
918
    if (lock_unlock(&finger->grav.mlock) != 0) error("Failed to unlock cell.");
919
920
921
922
923
924
925
926
927
928
929
930
931
932
  }

  /* 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)
933
      atomic_dec(&finger2->grav.mhold);
934
935

    /* Unlock this cell. */
936
    if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
937
938
939
940
941
942
943

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

944
945
946
947
948
949
950
951
952
953
954
/**
 * @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. */
955
  if (c->stars.hold || lock_trylock(&c->stars.lock) != 0) {
956
957
958
959
960
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
961
  if (c->stars.hold) {
962
963

    /* Unlock this cell. */
964
    if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
965
966
967
968
969
970
971
972
973
974
975

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

    /* Increment the hold. */
979
    atomic_inc(&finger->stars.hold);
980
981

    /* Unlock the cell. */
982
    if (lock_unlock(&finger->stars.lock) != 0) error("Failed to unlock cell.");
983
984
985
986
987
988
989
990
991
992
993
994
995
996
  }

  /* 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)
997
      atomic_dec(&finger2->stars.hold);
998
999

    /* Unlock this cell. */
1000
    if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
1001
1002
1003
1004
1005
1006
1007

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

1008
/**
1009
 * @brief Unlock a cell's parents for access to #part array.
1010
1011
1012
 *
 * @param c The #cell.
 */
1013
1014
1015
1016
1017
void cell_unlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
1021
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
1022
    atomic_dec(&finger->hydro.hold);
1023
1024
1025
1026

  TIMER_TOC(timer_locktree);
}

1027
1028
1029
1030
1031
/**
 * @brief Unlock a cell's parents for access to #gpart array.
 *
 * @param c The #cell.
 */
1032
1033
1034
1035
1036
void cell_gunlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
1040
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
1041
    atomic_dec(&finger->grav.phold);
1042
1043
1044
1045

  TIMER_TOC(timer_locktree);
}

1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
/**
 * @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. */
1056
  if (lock_unlock(&c->grav.mlock) != 0) error("Failed to unlock cell.");
1057
1058
1059

  /* Climb up the tree and unhold the parents. */
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
1060
    atomic_dec(&finger->grav.mhold);
1061
1062
1063
1064

  TIMER_TOC(timer_locktree);
}

1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
/**
 * @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. */
1075
  if (lock_unlock(&c->stars.lock) != 0) error("Failed to unlock cell.");
1076
1077
1078

  /* Climb up the tree and unhold the parents. */
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
1079
    atomic_dec(&finger->stars.hold);
1080
1081
1082
1083

  TIMER_TOC(timer_locktree);
}

1084
1085
1086
1087
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
1088
 * @param parts_offset Offset of the cell parts array relative to the
1089
 *        space's parts array, i.e. c->hydro.parts - s->parts.
1090
 * @param sparts_offset Offset of the cell sparts array relative to the
1091
1092
 *        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)
1093
 * entries, used for sorting indices.
1094
1095
1096
 * @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)
1097
 * entries, used for sorting indices for the gparts.
1098
 */
1099
1100
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
                struct cell_buff *buff, struct cell_buff *sbuff,
1101
                struct cell_buff *gbuff) {
1102

1103
1104
  const int count = c->hydro.count, gcount = c->grav.count,
            scount = c->stars.count;
1105
1106
  struct part *parts = c->hydro.parts;
  struct xpart *xparts = c->hydro.xparts;
1107
1108
  struct gpart *gparts = c->grav.parts;
  struct spart *sparts = c->stars.parts;
1109
1110
1111
1112
1113
1114
  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];

1115
1116
1117
#ifdef SWIFT_DEBUG_CHECKS
  /* Check that the buffs are OK. */
  for (int k = 0; k < count; k++) {
1118
    if (buff[k].x[0] != parts[k].x[0] || buff[k].x[1] != parts[k].x[1] ||
1119
        buff[k].x[2] != parts[k].x[2])
1120
1121
      error("Inconsistent buff contents.");
  }
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
  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.");
  }
1132
#endif /* SWIFT_DEBUG_CHECKS */
1133
1134
1135

  /* Fill the buffer with the indices. */
  for (int k = 0; k < count; k++) {
1136
1137
    const int bid = (buff[k].x[0] >= pivot[0]) * 4 +
                    (buff[k].x[1] >= pivot[1]) * 2 + (buff[k].x[2] >= pivot[2]);
1138
    bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
1139
    buff[k].ind = bid;
1140
  }
1141

1142
1143
1144
1145
1146
  /* 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;
1147
1148
  }

1149
1150
1151
1152
  /* 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
1153
      int bid = buff[k].ind;
1154