cell.c 73.7 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 "drift.h"
53
#include "error.h"
54
#include "gravity.h"
55
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
56
#include "hydro_properties.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
57
#include "memswap.h"
58
#include "minmax.h"
59
#include "scheduler.h"
60
61
#include "space.h"
#include "timers.h"
62

63
64
65
/* Global variables. */
int cell_next_tag = 0;

66
67
68
69
70
/**
 * @brief Get the size of the cell subtree.
 *
 * @param c The #cell.
 */
71
int cell_getsize(struct cell *c) {
72

Pedro Gonnet's avatar
Pedro Gonnet committed
73
74
  /* Number of cells in this subtree. */
  int count = 1;
75

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

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

/**
86
87
88
89
90
91
92
93
 * @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.
 *
 * @return The number of cells created.
 */
94
95
int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {

96
97
#ifdef WITH_MPI

98
99
  /* Unpack the current pcell. */
  c->h_max = pc->h_max;
100
101
  c->ti_end_min = pc->ti_end_min;
  c->ti_end_max = pc->ti_end_max;
102
103
  c->ti_old_part = pc->ti_old_part;
  c->ti_old_gpart = pc->ti_old_gpart;
104
  c->count = pc->count;
105
  c->gcount = pc->gcount;
106
  c->scount = pc->scount;
107
  c->tag = pc->tag;
Matthieu Schaller's avatar
Matthieu Schaller committed
108

109
110
  /* Number of new cells created. */
  int count = 1;
111
112

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
113
  for (int k = 0; k < 8; k++)
114
    if (pc->progeny[k] >= 0) {
115
116
      struct cell *temp;
      space_getcells(s, 1, &temp);
117
      temp->count = 0;
118
      temp->gcount = 0;
119
      temp->scount = 0;
120
121
122
      temp->loc[0] = c->loc[0];
      temp->loc[1] = c->loc[1];
      temp->loc[2] = c->loc[2];
123
124
125
      temp->width[0] = c->width[0] / 2;
      temp->width[1] = c->width[1] / 2;
      temp->width[2] = c->width[2] / 2;
126
      temp->dmin = c->dmin / 2;
127
128
129
      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];
130
131
      temp->depth = c->depth + 1;
      temp->split = 0;
132
133
      temp->dx_max_part = 0.f;
      temp->dx_max_gpart = 0.f;
134
      temp->dx_max_sort = 0.f;
135
136
137
138
139
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
      count += cell_unpack(&pc[pc->progeny[k]], temp, s);
140
141
    }

142
  /* Return the total number of unpacked cells. */
143
  c->pcell_size = count;
144
  return count;
145
146
147
148
149

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

152
/**
153
 * @brief Link the cells recursively to the given #part array.
154
155
156
157
158
159
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */
160
int cell_link_parts(struct cell *c, struct part *parts) {
161

162
163
164
  c->parts = parts;

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
165
166
167
168
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
169
        offset += cell_link_parts(c->progeny[k], &parts[offset]);
Pedro Gonnet's avatar
Pedro Gonnet committed
170
171
    }
  }
172

173
  /* Return the total number of linked particles. */
174
175
  return c->count;
}
176

177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
/**
 * @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) {

  c->gparts = gparts;

  /* 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. */
  return c->gcount;
}

202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
/**
 * @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) {

  c->sparts = sparts;

  /* 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. */
  return c->scount;
}

227
228
229
230
231
232
233
234
235
/**
 * @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.
 *
 * @return The number of packed cells.
 */
236
237
int cell_pack(struct cell *c, struct pcell *pc) {

238
239
#ifdef WITH_MPI

240
241
  /* Start by packing the data of the current cell. */
  pc->h_max = c->h_max;
242
243
  pc->ti_end_min = c->ti_end_min;
  pc->ti_end_max = c->ti_end_max;
244
245
  pc->ti_old_part = c->ti_old_part;
  pc->ti_old_gpart = c->ti_old_gpart;
246
  pc->count = c->count;
247
  pc->gcount = c->gcount;
248
  pc->scount = c->scount;
249
250
251
  c->tag = pc->tag = atomic_inc(&cell_next_tag) % cell_max_tag;

  /* Fill in the progeny, depth-first recursion. */
Pedro Gonnet's avatar
Pedro Gonnet committed
252
253
  int count = 1;
  for (int k = 0; k < 8; k++)
254
255
256
257
258
259
260
    if (c->progeny[k] != NULL) {
      pc->progeny[k] = count;
      count += cell_pack(c->progeny[k], &pc[count]);
    } else
      pc->progeny[k] = -1;

  /* Return the number of packed cells used. */
261
262
  c->pcell_size = count;
  return count;
263
264
265
266
267

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

270
271
272
273
274
275
276
277
/**
 * @brief Pack the time information of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
 * @param ti_ends (output) The time information we pack into
 *
 * @return The number of packed cells.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
278
int cell_pack_ti_ends(struct cell *c, integertime_t *ti_ends) {
279

280
281
#ifdef WITH_MPI

282
283
  /* Pack this cell's data. */
  ti_ends[0] = c->ti_end_min;
284

285
286
287
288
289
290
291
292
293
  /* 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_ti_ends(c->progeny[k], &ti_ends[count]);
    }

  /* Return the number of packed values. */
  return count;
294
295
296
297
298

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

301
302
303
304
/**
 * @brief Unpack the time information of a given cell and its sub-cells.
 *
 * @param c The #cell
305
 * @param ti_ends The time information to unpack
306
307
308
 *
 * @return The number of cells created.
 */
309
int cell_unpack_ti_ends(struct cell *c, integertime_t *ti_ends) {
310

311
312
#ifdef WITH_MPI

313
314
  /* Unpack this cell's data. */
  c->ti_end_min = ti_ends[0];
315

316
317
318
319
320
321
322
323
  /* 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_ti_ends(c->progeny[k], &ti_ends[count]);
    }

  /* Return the number of packed values. */
324
  return count;
325
326
327
328
329

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

332
/**
333
 * @brief Lock a cell for access to its array of #part and hold its parents.
334
335
 *
 * @param c The #cell.
336
 * @return 0 on success, 1 on failure
337
 */
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
int cell_locktree(struct cell *c) {

  TIMER_TIC

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

  /* Did somebody hold this cell in the meantime? */
  if (c->hold) {

    /* Unlock this cell. */
    if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

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

  /* Climb up the tree and lock/hold/unlock. */
Pedro Gonnet's avatar
Pedro Gonnet committed
360
  struct cell *finger;
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

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

    /* Increment the hold. */
    atomic_inc(&finger->hold);

    /* Unlock the cell. */
    if (lock_unlock(&finger->lock) != 0) error("Failed to unlock cell.");
  }

  /* 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
383
384
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
385
      atomic_dec(&finger2->hold);
386
387
388
389
390
391
392
393
394
395

    /* Unlock this cell. */
    if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

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

396
397
398
399
400
401
/**
 * @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
 */
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
int cell_glocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
  if (c->ghold || lock_trylock(&c->glock) != 0) {
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
  if (c->ghold) {

    /* Unlock this cell. */
    if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

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

  /* Climb up the tree and lock/hold/unlock. */
Pedro Gonnet's avatar
Pedro Gonnet committed
424
  struct cell *finger;
425
426
427
428
429
430
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

    /* Lock this cell. */
    if (lock_trylock(&finger->glock) != 0) break;

    /* Increment the hold. */
431
    atomic_inc(&finger->ghold);
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446

    /* Unlock the cell. */
    if (lock_unlock(&finger->glock) != 0) error("Failed to unlock cell.");
  }

  /* 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
447
448
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
449
      atomic_dec(&finger2->ghold);
450
451
452
453
454
455
456
457
458

    /* Unlock this cell. */
    if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

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

460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
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
512
513
514
515
516
517
518
519
520
521
522
523
/**
 * @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. */
  if (c->mhold || lock_trylock(&c->mlock) != 0) {
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
  if (c->mhold) {

    /* Unlock this cell. */
    if (lock_unlock(&c->mlock) != 0) error("Failed to unlock cell.");

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

    /* Increment the hold. */
    atomic_inc(&finger->mhold);

    /* Unlock the cell. */
    if (lock_unlock(&finger->mlock) != 0) error("Failed to unlock cell.");
  }

  /* 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)
      atomic_dec(&finger2->mhold);

    /* Unlock this cell. */
    if (lock_unlock(&c->mlock) != 0) error("Failed to unlock cell.");

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

524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
/**
 * @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. */
  if (c->shold || lock_trylock(&c->slock) != 0) {
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
  if (c->shold) {

    /* Unlock this cell. */
    if (lock_unlock(&c->slock) != 0) error("Failed to unlock cell.");

    /* 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. */
    if (lock_trylock(&finger->slock) != 0) break;

    /* Increment the hold. */
    atomic_inc(&finger->shold);

    /* Unlock the cell. */
    if (lock_unlock(&finger->slock) != 0) error("Failed to unlock cell.");
  }

  /* 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)
      atomic_dec(&finger2->shold);

    /* Unlock this cell. */
    if (lock_unlock(&c->slock) != 0) error("Failed to unlock cell.");

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

588
/**
589
 * @brief Unlock a cell's parents for access to #part array.
590
591
592
 *
 * @param c The #cell.
 */
593
594
595
596
597
598
599
600
void cell_unlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
601
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
602
    atomic_dec(&finger->hold);
603
604
605
606

  TIMER_TOC(timer_locktree);
}

607
608
609
610
611
/**
 * @brief Unlock a cell's parents for access to #gpart array.
 *
 * @param c The #cell.
 */
612
613
614
615
616
617
618
619
void cell_gunlocktree(struct cell *c) {

  TIMER_TIC

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

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
620
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
621
    atomic_dec(&finger->ghold);
622
623
624
625

  TIMER_TOC(timer_locktree);
}

626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
/**
 * @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. */
  if (lock_unlock(&c->mlock) != 0) error("Failed to unlock cell.");

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

  TIMER_TOC(timer_locktree);
}

645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
/**
 * @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. */
  if (lock_unlock(&c->slock) != 0) error("Failed to unlock cell.");

  /* Climb up the tree and unhold the parents. */
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
    atomic_dec(&finger->shold);

  TIMER_TOC(timer_locktree);
}

664
665
666
667
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
668
669
 * @param parts_offset Offset of the cell parts array relative to the
 *        space's parts array, i.e. c->parts - s->parts.
670
671
 * @param sparts_offset Offset of the cell sparts array relative to the
 *        space's sparts array, i.e. c->sparts - s->sparts.
672
673
 * @param buff A buffer with at least max(c->count, c->gcount) entries,
 *        used for sorting indices.
674
675
 * @param sbuff A buffer with at least max(c->scount, c->gcount) entries,
 *        used for sorting indices for the sparts.
Peter W. Draper's avatar
Peter W. Draper committed
676
677
 * @param gbuff A buffer with at least max(c->count, c->gcount) entries,
 *        used for sorting indices for the gparts.
678
 */
679
680
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
                struct cell_buff *buff, struct cell_buff *sbuff,
681
                struct cell_buff *gbuff) {
682

683
  const int count = c->count, gcount = c->gcount, scount = c->scount;
684
685
686
  struct part *parts = c->parts;
  struct xpart *xparts = c->xparts;
  struct gpart *gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
687
  struct spart *sparts = c->sparts;
688
689
690
691
692
693
  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];

694
695
696
#ifdef SWIFT_DEBUG_CHECKS
  /* Check that the buffs are OK. */
  for (int k = 0; k < count; k++) {
697
    if (buff[k].x[0] != parts[k].x[0] || buff[k].x[1] != parts[k].x[1] ||
698
        buff[k].x[2] != parts[k].x[2])
699
700
      error("Inconsistent buff contents.");
  }
701
702
703
704
705
706
707
708
709
710
  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.");
  }
711
#endif /* SWIFT_DEBUG_CHECKS */
712
713
714

  /* Fill the buffer with the indices. */
  for (int k = 0; k < count; k++) {
715
716
    const int bid = (buff[k].x[0] > pivot[0]) * 4 +
                    (buff[k].x[1] > pivot[1]) * 2 + (buff[k].x[2] > pivot[2]);
717
    bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
718
    buff[k].ind = bid;
719
  }
720

721
722
723
724
725
  /* 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;
726
727
  }

728
729
730
731
  /* 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
732
      int bid = buff[k].ind;
733
734
735
      if (bid != bucket) {
        struct part part = parts[k];
        struct xpart xpart = xparts[k];
736
        struct cell_buff temp_buff = buff[k];
737
738
        while (bid != bucket) {
          int j = bucket_offset[bid] + bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
739
          while (buff[j].ind == bid) {
740
741
742
            j++;
            bucket_count[bid]++;
          }
Pedro Gonnet's avatar
Pedro Gonnet committed
743
744
          memswap(&parts[j], &part, sizeof(struct part));
          memswap(&xparts[j], &xpart, sizeof(struct xpart));
745
746
          memswap(&buff[j], &temp_buff, sizeof(struct cell_buff));
          bid = temp_buff.ind;
747
748
749
        }
        parts[k] = part;
        xparts[k] = xpart;
750
        buff[k] = temp_buff;
751
      }
752
      bucket_count[bid]++;
753
754
755
756
    }
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
757
  for (int k = 0; k < 8; k++) {
758
759
760
    c->progeny[k]->count = bucket_count[k];
    c->progeny[k]->parts = &c->parts[bucket_offset[k]];
    c->progeny[k]->xparts = &c->xparts[bucket_offset[k]];
761
762
763
  }

  /* Re-link the gparts. */
764
765
  if (count > 0 && gcount > 0)
    part_relink_gparts_to_parts(parts, count, parts_offset);
766

767
#ifdef SWIFT_DEBUG_CHECKS
768
769
770
771
772
773
774
775
  /* 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);
  }

776
  /* Verify that _all_ the parts have been assigned to a cell. */
777
778
779
780
781
782
783
784
  for (int k = 1; k < 8; k++)
    if (&c->progeny[k - 1]->parts[c->progeny[k - 1]->count] !=
        c->progeny[k]->parts)
      error("Particle sorting failed (internal consistency).");
  if (c->progeny[0]->parts != c->parts)
    error("Particle sorting failed (left edge).");
  if (&c->progeny[7]->parts[c->progeny[7]->count] != &c->parts[count])
    error("Particle sorting failed (right edge).");
785
786

  /* Verify a few sub-cells. */
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
  for (int k = 0; k < c->progeny[0]->count; k++)
    if (c->progeny[0]->parts[k].x[0] > pivot[0] ||
        c->progeny[0]->parts[k].x[1] > pivot[1] ||
        c->progeny[0]->parts[k].x[2] > pivot[2])
      error("Sorting failed (progeny=0).");
  for (int k = 0; k < c->progeny[1]->count; k++)
    if (c->progeny[1]->parts[k].x[0] > pivot[0] ||
        c->progeny[1]->parts[k].x[1] > pivot[1] ||
        c->progeny[1]->parts[k].x[2] <= pivot[2])
      error("Sorting failed (progeny=1).");
  for (int k = 0; k < c->progeny[2]->count; k++)
    if (c->progeny[2]->parts[k].x[0] > pivot[0] ||
        c->progeny[2]->parts[k].x[1] <= pivot[1] ||
        c->progeny[2]->parts[k].x[2] > pivot[2])
      error("Sorting failed (progeny=2).");
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
  for (int k = 0; k < c->progeny[3]->count; k++)
    if (c->progeny[3]->parts[k].x[0] > pivot[0] ||
        c->progeny[3]->parts[k].x[1] <= pivot[1] ||
        c->progeny[3]->parts[k].x[2] <= pivot[2])
      error("Sorting failed (progeny=3).");
  for (int k = 0; k < c->progeny[4]->count; k++)
    if (c->progeny[4]->parts[k].x[0] <= pivot[0] ||
        c->progeny[4]->parts[k].x[1] > pivot[1] ||
        c->progeny[4]->parts[k].x[2] > pivot[2])
      error("Sorting failed (progeny=4).");
  for (int k = 0; k < c->progeny[5]->count; k++)
    if (c->progeny[5]->parts[k].x[0] <= pivot[0] ||
        c->progeny[5]->parts[k].x[1] > pivot[1] ||
        c->progeny[5]->parts[k].x[2] <= pivot[2])
      error("Sorting failed (progeny=5).");
  for (int k = 0; k < c->progeny[6]->count; k++)
    if (c->progeny[6]->parts[k].x[0] <= pivot[0] ||
        c->progeny[6]->parts[k].x[1] <= pivot[1] ||
        c->progeny[6]->parts[k].x[2] > pivot[2])
      error("Sorting failed (progeny=6).");
  for (int k = 0; k < c->progeny[7]->count; k++)
    if (c->progeny[7]->parts[k].x[0] <= pivot[0] ||
        c->progeny[7]->parts[k].x[1] <= pivot[1] ||
        c->progeny[7]->parts[k].x[2] <= pivot[2])
      error("Sorting failed (progeny=7).");
827
#endif
828

829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
  /* 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));
          bid = temp_buff.ind;
        }
        sparts[k] = spart;
        sbuff[k] = temp_buff;
      }
      bucket_count[bid]++;
    }
  }

  /* Store the counts and offsets. */
  for (int k = 0; k < 8; k++) {
    c->progeny[k]->scount = bucket_count[k];
    c->progeny[k]->sparts = &c->sparts[bucket_offset[k]];
  }

  /* Re-link the gparts. */
  if (scount > 0 && gcount > 0)
880
    part_relink_gparts_to_sparts(sparts, scount, sparts_offset);
881
882

  /* Finally, do the same song and dance for the gparts. */
883
884
885
886
  for (int k = 0; k < 8; k++) bucket_count[k] = 0;

  /* Fill the buffer with the indices. */
  for (int k = 0; k < gcount; k++) {
887
888
    const int bid = (gbuff[k].x[0] > pivot[0]) * 4 +
                    (gbuff[k].x[1] > pivot[1]) * 2 + (gbuff[k].x[2] > pivot[2]);
889
    bucket_count[bid]++;
890
    gbuff[k].ind = bid;
891
  }
892
893
894
895
896
897

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

900
901
902
903
  /* 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++) {
904
      int bid = gbuff[k].ind;
905
906
      if (bid != bucket) {
        struct gpart gpart = gparts[k];
907
        struct cell_buff temp_buff = gbuff[k];
908
909
        while (bid != bucket) {
          int j = bucket_offset[bid] + bucket_count[bid]++;
910
          while (gbuff[j].ind == bid) {
911
912
913
            j++;
            bucket_count[bid]++;
          }
Pedro Gonnet's avatar
Pedro Gonnet committed
914
          memswap(&gparts[j], &gpart, sizeof(struct gpart));
915
916
          memswap(&gbuff[j], &temp_buff, sizeof(struct cell_buff));
          bid = temp_buff.ind;
917
918
        }
        gparts[k] = gpart;
919
        gbuff[k] = temp_buff;
920
      }
921
      bucket_count[bid]++;
922
923
924
925
    }
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
926
  for (int k = 0; k < 8; k++) {
927
928
    c->progeny[k]->gcount = bucket_count[k];
    c->progeny[k]->gparts = &c->gparts[bucket_offset[k]];
929
930
931
  }

  /* Re-link the parts. */
932
  if (count > 0 && gcount > 0)
933
    part_relink_parts_to_gparts(gparts, gcount, parts - parts_offset);
934
935
936
937

  /* Re-link the sparts. */
  if (scount > 0 && gcount > 0)
    part_relink_sparts_to_gparts(gparts, gcount, sparts - sparts_offset);
938
}
939

940
941
942
943
/**
 * @brief Sanitizes the smoothing length values of cells by setting large
 * outliers to more sensible values.
 *
944
945
 * Each cell with <1000 part will be processed. We limit h to be the size of
 * the cell and replace 0s with a good estimate.
946
947
 *
 * @param c The cell.
948
 * @param treated Has the cell already been sanitized at this level ?
949
 */
950
void cell_sanitize(struct cell *c, int treated) {
951
952
953

  const int count = c->count;
  struct part *parts = c->parts;
954
  float h_max = 0.f;
955

956
957
  /* Treat cells will <1000 particles */
  if (count < 1000 && !treated) {
958

959
960
    /* Get an upper bound on h */
    const float upper_h_max = c->dmin / (1.2f * kernel_gamma);
961

962
963
964
965
966
967
    /* Apply it */
    for (int i = 0; i < count; ++i) {
      if (parts[i].h == 0.f || parts[i].h > upper_h_max)
        parts[i].h = upper_h_max;
    }
  }
968

969
970
  /* Recurse and gather the new h_max values */
  if (c->split) {
971

972
973
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
974

975
976
        /* Recurse */
        cell_sanitize(c->progeny[k], (count < 1000));
977

978
979
980
981
        /* And collect */
        h_max = max(h_max, c->progeny[k]->h_max);
      }
    }
982
983
  } else {

984
985
    /* Get the new value of h_max */
    for (int i = 0; i < count; ++i) h_max = max(h_max, parts[i].h);
986
  }
987
988
989

  /* Record the change */
  c->h_max = h_max;
990
991
}

992
/**
993
 * @brief Converts hydro quantities to a valid state after the initial density
994
 * calculation
995
996
997
998
999
1000
1001
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_convert_hydro(struct cell *c, void *data) {

  struct part *p = c->parts;
1002
  struct xpart *xp = c->xparts;
1003
1004

  for (int i = 0; i < c->count; ++i) {
1005
    hydro_convert_quantities(&p[i], &xp[i]);
1006
1007
1008
  }
}

Matthieu Schaller's avatar
Matthieu Schaller committed
1009
1010
1011
1012
1013
1014
/**
 * @brief Cleans the links in a given cell.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
1015
void cell_clean_links(struct cell *c, void *data) {
Matthieu Schaller's avatar
Matthieu Schaller committed
1016
  c->density = NULL;
1017
  c->gradient = NULL;
Matthieu Schaller's avatar
Matthieu Schaller committed
1018
  c->force = NULL;
1019
  c->grav = NULL;
Matthieu Schaller's avatar
Matthieu Schaller committed
1020
}
1021

1022
/**
1023
 * @brief Checks that the #part in a cell are at the
1024
 * current point in time
1025
1026
1027
1028
1029
1030
 *
 * Calls error() if the cell is not at the current time.
 *
 * @param c Cell to act upon
 * @param data The current time on the integer time-line
 */
1031
void cell_check_part_drift_point(struct cell *c, void *data) {
1032

1033
1034
#ifdef SWIFT_DEBUG_CHECKS

1035
  const integertime_t ti_drift = *(integertime_t *)data;
1036

1037
  /* Only check local cells */
1038
  if (c->nodeID != engine_rank) return;
1039

1040
1041
1042
  if (c->ti_old_part != ti_drift)
    error("Cell in an incorrect time-zone! c->ti_old_part=%lld ti_drift=%lld",
          c->ti_old_part, ti_drift);
1043

1044
1045
  for (int i = 0; i < c->count; ++i)
    if (c->parts[i].ti_drift != ti_drift)
1046
      error("part in an incorrect time-zone! p->ti_drift=%lld ti_drift=%lld",
1047
            c->parts[i].ti_drift, ti_drift);
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

/**
 * @brief Checks that the #gpart and #spart in a cell are at the
 * current point in time
 *
 * Calls error() if the cell is not at the current time.
 *
 * @param c Cell to act upon
 * @param data The current time on the integer time-line
 */
void cell_check_gpart_drift_point(struct cell *c, void *data) {

#ifdef SWIFT_DEBUG_CHECKS

  const integertime_t ti_drift = *(integertime_t *)data;

  /* Only check local cells */
  if (c->nodeID != engine_rank) return;

  if (c->ti_old_gpart != ti_drift)
    error("Cell in an incorrect time-zone! c->ti_old_gpart=%lld ti_drift=%lld",
          c->ti_old_gpart, ti_drift);
1074

1075
1076
  for (int i = 0; i < c->gcount; ++i)
    if (c->gparts[i].ti_drift != ti_drift)
1077
      error("g-part in an incorrect time-zone! gp->ti_drift=%lld ti_drift=%lld",
1078
            c->gparts[i].ti_drift, ti_drift);
1079

1080
1081
1082
1083
  for (int i = 0; i < c->scount; ++i)
    if (c->sparts[i].ti_drift != ti_drift)
      error("s-part in an incorrect time-zone! sp->ti_drift=%lld ti_drift=%lld",
            c->sparts[i].ti_drift, ti_drift);
1084
1085
1086
#else
  error("Calling debugging code without debugging flag activated.");
#endif
1087
1088
}

1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
/**
 * @brief Checks that the multipole of a cell is at the current point in time
 *
 * Calls error() if the cell is not at the current time.
 *
 * @param c Cell to act upon
 * @param data The current time on the integer time-line
 */
void cell_check_multipole_drift_point(struct cell *c, void *data) {

#ifdef SWIFT_DEBUG_CHECKS

  const integertime_t ti_drift = *(integertime_t *)data;

  /* Only check local cells */
  if (c->nodeID != engine_rank) return;

  if (c->ti_old_multipole != ti_drift)
    error(
        "Cell multipole in an incorrect time-zone! c->ti_old_multipole=%lld "
1109
1110
        "ti_drift=%lld (depth=%d)",
        c->ti_old_multipole, ti_drift, c->depth);
1111
1112
1113
1114
1115
1116

#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
/**
 * @brief Resets all the individual cell task counters to 0.
 *
 * Should only be used for debugging purposes.
 *
 * @param c The #cell to reset.
 */
void cell_reset_task_counters(struct cell *c) {

#ifdef SWIFT_DEBUG_CHECKS
  for (int t = 0; t < task_type_count; ++t) c->tasks_executed[t] = 0;
  for (int t = 0; t < task_subtype_count; ++t) c->subtasks_executed[t] = 0;
#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

1134
1135
1136
1137
/**
 * @brief Recursively construct all the multipoles in a cell hierarchy.
 *
 * @param c The #cell.
Matthieu Schaller's avatar
Matthieu Schaller committed
1138
 * @param ti_current The current integer time.
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
 */
void cell_make_multipoles(struct cell *c, integertime_t ti_current) {

  /* Reset everything */
  gravity_reset(c->multipole);

  if (c->split) {

    /* Compute CoM of all progenies */
    double CoM[3] = {0., 0., 0.};
    double mass = 0.;

    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
        const struct gravity_tensors *m = c->progeny[k]->multipole;
        CoM[0] += m->CoM[0] * m->m_pole.M_000;
        CoM[1] += m->CoM[1] * m->m_pole.M_000;
        CoM[2] += m->CoM[2] * m->m_pole.M_000;
        mass += m->m_pole.M_000;
      }
    }
    c->multipole->CoM[0] = CoM[0] / mass;
    c->multipole->CoM[1] = CoM[1] / mass;
    c->multipole->CoM[2] = CoM[2] / mass;

    /* Now shift progeny multipoles and add them up */
    struct multipole temp;
    double r_max = 0.;
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
        const struct cell *cp = c->progeny[k];
        const struct multipole *m = &cp->multipole->m_pole;

        /* Contribution to multipole */
        gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM);
        gravity_multipole_add(&c->multipole->m_pole, &temp);

        /* Upper limit of max CoM<->gpart distance */
        const double dx = c->multipole->CoM[0] - cp->multipole->CoM[0];
        const double dy = c->multipole->CoM[1] - cp->multipole->CoM[1];
        const double dz = c->multipole->CoM[2] - cp->multipole->CoM[2];
        const double r2 = dx * dx + dy * dy + dz * dz;
        r_max = max(r_max, cp->multipole->r_max + sqrt(r2));
      }
    }
    /* Alternative upper limit of max CoM<->gpart distance */
    const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
                          ? c->multipole->CoM[0] - c->loc[0]
                          : c->loc[0] + c->width[0] - c->multipole->CoM[0];
    const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
                          ? c->multipole->CoM[1] - c->loc[1]
                          : c->loc[1] + c->width[1] - c->multipole->CoM[1];
    const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
                          ? c->multipole->CoM[2] - c->loc[2]
                          : c->loc[2] + c->width[2] - c->multipole->CoM[2];

    /* Take minimum of both limits */
    c->multipole->r_max = min(r_max, sqrt(dx * dx + dy * dy + dz * dz));

  } else {

    if (c->gcount > 0) {
      gravity_P2M(c->multipole, c->gparts, c->gcount);
      const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
                            ? c->multipole->CoM[0] - c->loc[0]
                            : c->loc[0] + c->width[0] - c->multipole->CoM[0];
      const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
                            ? c->multipole->CoM[1] - c->loc[1]
                            : c->loc[1] + c->width[1] - c->multipole->CoM[1];
      const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
                            ? c->multipole->CoM[2] - c->loc[2]
                            : c->loc[2] + c->width[2] - c->multipole->CoM[2];
      c->multipole->r_max = sqrt(dx * dx + dy * dy + dz * dz);
    } else {
      gravity_multipole_init(&c->multipole->m_pole);
      c->multipole->CoM[0] = c->loc[0] + c->width[0] / 2.;
      c->multipole->CoM[1] = c->loc[1] + c->width[1] / 2.;
      c->multipole->CoM[2] = c->loc[2] + c->width[2] / 2.;
      c->multipole->r_max = 0.;
    }
  }

  c->ti_old_multipole = ti_current;
}

1224
1225
1226
1227
1228
1229
1230
1231
1232
/**
 * @brief Computes the multi-pole brutally and compare to the
 * recursively computed one.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_check_multipole(struct cell *c, void *data) {

1233
#ifdef SWIFT_DEBUG_CHECKS
1234
  struct gravity_tensors ma;
1235
  const double tolerance = 1e-3; /* Relative */
1236

1237
1238
  return;

1239
1240
1241
1242
  /* First recurse */
  if (c->split)
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) cell_check_multipole(c->progeny[k], NULL);
1243
1244
1245
1246

  if (c->gcount > 0) {

    /* Brute-force calculation */
1247
    gravity_P2M(&ma, c->gparts, c->gcount);
1248
1249

    /* Now  compare the multipole expansion */
1250
    if (!gravity_multipole_equal(&ma, c->multipole, tolerance)) {
Matthieu Schaller's avatar
Matthieu Schaller committed
1251
1252
      message("Multipoles are not equal at depth=%d! tol=%f", c->depth,
              tolerance);
1253
      message("Correct answer:");
1254
      gravity_multipole_print(&ma.m_pole);
1255
      message("Recursive multipole:");
1256
      gravity_multipole_print(&c->multipole->m_pole);
1257
1258
      error("Aborting");
    }
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268

    /* Check that the upper limit of r_max is good enough */
    if (!(c->multipole->r_max >= ma.r_max)) {
      error("Upper-limit r_max=%e too small. Should be >=%e.",
            c->multipole->r_max, ma.r_max);
    } else if (c->multipole->r_max * c->multipole->r_max >
               3. * c->width[0] * c->width[0]) {
      error("r_max=%e larger than cell diagonal %e.", c->multipole->r_max,
            sqrt(3. * c->width[0] * c->width[0]));
    }
Matthieu Schaller's avatar