cell.c 84.2 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 63
#include "space.h"
#include "timers.h"
64

65 66 67
/* Global variables. */
int cell_next_tag = 0;

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

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

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

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

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

97 98 99
  c->parts = parts;

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
100 101 102 103
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
104
        offset += cell_link_parts(c->progeny[k], &parts[offset]);
Pedro Gonnet's avatar
Pedro Gonnet committed
105 106
    }
  }
107

108
  /* Return the total number of linked particles. */
109 110
  return c->count;
}
111

112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
/**
 * @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;
}

137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
/**
 * @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;
}

162 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.
 *
 * @return The number of packed cells.
 */
171
int cell_pack(struct cell *restrict c, struct pcell *restrict pc) {
172

173 174
#ifdef WITH_MPI

175 176
  /* Start by packing the data of the current cell. */
  pc->h_max = c->h_max;
177 178 179 180
  pc->ti_hydro_end_min = c->ti_hydro_end_min;
  pc->ti_hydro_end_max = c->ti_hydro_end_max;
  pc->ti_gravity_end_min = c->ti_gravity_end_min;
  pc->ti_gravity_end_max = c->ti_gravity_end_max;
181 182
  pc->ti_old_part = c->ti_old_part;
  pc->ti_old_gpart = c->ti_old_gpart;
183
  pc->ti_old_multipole = c->ti_old_multipole;
184
  pc->count = c->count;
185
  pc->gcount = c->gcount;
186
  pc->scount = c->scount;
187
  c->tag = pc->tag = atomic_inc(&cell_next_tag) % cell_max_tag;
188 189 190
#ifdef SWIFT_DEBUG_CHECKS
  pc->cellID = c->cellID;
#endif
191 192

  /* Fill in the progeny, depth-first recursion. */
Pedro Gonnet's avatar
Pedro Gonnet committed
193 194
  int count = 1;
  for (int k = 0; k < 8; k++)
195 196 197 198 199 200 201
    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. */
202 203
  c->pcell_size = count;
  return count;
204 205 206 207 208

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

211 212 213 214 215 216 217 218 219
/**
 * @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.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
220 221
int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
                struct space *restrict s) {
222 223 224 225 226

#ifdef WITH_MPI

  /* Unpack the current pcell. */
  c->h_max = pc->h_max;
227 228 229 230
  c->ti_hydro_end_min = pc->ti_hydro_end_min;
  c->ti_hydro_end_max = pc->ti_hydro_end_max;
  c->ti_gravity_end_min = pc->ti_gravity_end_min;
  c->ti_gravity_end_max = pc->ti_gravity_end_max;
231 232
  c->ti_old_part = pc->ti_old_part;
  c->ti_old_gpart = pc->ti_old_gpart;
233
  c->ti_old_multipole = pc->ti_old_multipole;
234 235 236 237
  c->count = pc->count;
  c->gcount = pc->gcount;
  c->scount = pc->scount;
  c->tag = pc->tag;
238 239 240
#ifdef SWIFT_DEBUG_CHECKS
  c->cellID = pc->cellID;
#endif
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284

  /* Number of new cells created. */
  int count = 1;

  /* Fill the progeny recursively, depth-first. */
  for (int k = 0; k < 8; k++)
    if (pc->progeny[k] >= 0) {
      struct cell *temp;
      space_getcells(s, 1, &temp);
      temp->count = 0;
      temp->gcount = 0;
      temp->scount = 0;
      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;
      temp->dx_max_part = 0.f;
      temp->dx_max_gpart = 0.f;
      temp->dx_max_sort = 0.f;
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
      count += cell_unpack(&pc[pc->progeny[k]], temp, s);
    }

  /* Return the total number of unpacked cells. */
  c->pcell_size = count;
  return count;

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

285 286 287 288
/**
 * @brief Pack the time information of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
289
 * @param pcells (output) The end-of-timestep information we pack into
290 291 292
 *
 * @return The number of packed cells.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
293 294
int cell_pack_end_step(struct cell *restrict c,
                       struct pcell_step *restrict pcells) {
295

296 297
#ifdef WITH_MPI

298
  /* Pack this cell's data. */
299
  pcells[0].ti_hydro_end_min = c->ti_hydro_end_min;
300
  pcells[0].ti_hydro_end_max = c->ti_hydro_end_max;
301
  pcells[0].ti_gravity_end_min = c->ti_gravity_end_min;
302
  pcells[0].ti_gravity_end_max = c->ti_gravity_end_max;
303 304
  pcells[0].dx_max_part = c->dx_max_part;
  pcells[0].dx_max_gpart = c->dx_max_gpart;
305

306 307 308 309
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
310
      count += cell_pack_end_step(c->progeny[k], &pcells[count]);
311 312 313 314
    }

  /* Return the number of packed values. */
  return count;
315 316 317 318 319

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

322 323 324 325
/**
 * @brief Unpack the time information of a given cell and its sub-cells.
 *
 * @param c The #cell
326
 * @param pcells The end-of-timestep information to unpack
327 328 329
 *
 * @return The number of cells created.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
330 331
int cell_unpack_end_step(struct cell *restrict c,
                         struct pcell_step *restrict pcells) {
332

333 334
#ifdef WITH_MPI

335
  /* Unpack this cell's data. */
336
  c->ti_hydro_end_min = pcells[0].ti_hydro_end_min;
337
  c->ti_hydro_end_max = pcells[0].ti_hydro_end_max;
338
  c->ti_gravity_end_min = pcells[0].ti_gravity_end_min;
339
  c->ti_gravity_end_max = pcells[0].ti_gravity_end_max;
340 341
  c->dx_max_part = pcells[0].dx_max_part;
  c->dx_max_gpart = pcells[0].dx_max_gpart;
342

343 344 345 346
  /* Fill in the progeny, depth-first recursion. */
  int count = 1;
  for (int k = 0; k < 8; k++)
    if (c->progeny[k] != NULL) {
347
      count += cell_unpack_end_step(c->progeny[k], &pcells[count]);
348 349 350
    }

  /* Return the number of packed values. */
351
  return count;
352 353 354 355 356

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

359
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
360 361
 * @brief Pack the multipole information of the given cell and all it's
 * sub-cells.
362 363 364 365 366 367 368
 *
 * @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
369
                         struct gravity_tensors *restrict pcells) {
370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400

#ifdef WITH_MPI

  /* Pack this cell's data. */
  pcells[0] = *c->multipole;

  /* 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
401
                           struct gravity_tensors *restrict pcells) {
402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423

#ifdef WITH_MPI

  /* Unpack this cell's data. */
  *c->multipole = pcells[0];

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

424
/**
425
 * @brief Lock a cell for access to its array of #part and hold its parents.
426 427
 *
 * @param c The #cell.
428
 * @return 0 on success, 1 on failure
429
 */
430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451
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
452
  struct cell *finger;
453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474
  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
475 476
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
477
      atomic_dec(&finger2->hold);
478 479 480 481 482 483 484 485 486 487

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

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

488 489 490 491 492 493
/**
 * @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
 */
494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515
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
516
  struct cell *finger;
517 518 519 520 521 522
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

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

    /* Increment the hold. */
523
    atomic_inc(&finger->ghold);
524 525 526 527 528 529 530 531 532 533 534 535 536 537 538

    /* 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
539 540
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
541
      atomic_dec(&finger2->ghold);
542 543 544 545 546 547 548 549 550

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

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }
}
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 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615
/**
 * @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;
  }
}

616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679
/**
 * @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;
  }
}

680
/**
681
 * @brief Unlock a cell's parents for access to #part array.
682 683 684
 *
 * @param c The #cell.
 */
685 686 687 688 689 690 691 692
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
693
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
694
    atomic_dec(&finger->hold);
695 696 697 698

  TIMER_TOC(timer_locktree);
}

699 700 701 702 703
/**
 * @brief Unlock a cell's parents for access to #gpart array.
 *
 * @param c The #cell.
 */
704 705 706 707 708 709 710 711
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
712
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
713
    atomic_dec(&finger->ghold);
714 715 716 717

  TIMER_TOC(timer_locktree);
}

718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736
/**
 * @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);
}

737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755
/**
 * @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);
}

756 757 758 759
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
760 761
 * @param parts_offset Offset of the cell parts array relative to the
 *        space's parts array, i.e. c->parts - s->parts.
762 763
 * @param sparts_offset Offset of the cell sparts array relative to the
 *        space's sparts array, i.e. c->sparts - s->sparts.
764 765
 * @param buff A buffer with at least max(c->count, c->gcount) entries,
 *        used for sorting indices.
766 767
 * @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
768 769
 * @param gbuff A buffer with at least max(c->count, c->gcount) entries,
 *        used for sorting indices for the gparts.
770
 */
771 772
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
                struct cell_buff *buff, struct cell_buff *sbuff,
773
                struct cell_buff *gbuff) {
774

775
  const int count = c->count, gcount = c->gcount, scount = c->scount;
776 777 778
  struct part *parts = c->parts;
  struct xpart *xparts = c->xparts;
  struct gpart *gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
779
  struct spart *sparts = c->sparts;
780 781 782 783 784 785
  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];

786 787 788
#ifdef SWIFT_DEBUG_CHECKS
  /* Check that the buffs are OK. */
  for (int k = 0; k < count; k++) {
789
    if (buff[k].x[0] != parts[k].x[0] || buff[k].x[1] != parts[k].x[1] ||
790
        buff[k].x[2] != parts[k].x[2])
791 792
      error("Inconsistent buff contents.");
  }
793 794 795 796 797 798 799 800 801 802
  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.");
  }
803
#endif /* SWIFT_DEBUG_CHECKS */
804 805 806

  /* Fill the buffer with the indices. */
  for (int k = 0; k < count; k++) {
807 808
    const int bid = (buff[k].x[0] >= pivot[0]) * 4 +
                    (buff[k].x[1] >= pivot[1]) * 2 + (buff[k].x[2] >= pivot[2]);
809
    bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
810
    buff[k].ind = bid;
811
  }
812

813 814 815 816 817
  /* 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;
818 819
  }

820 821 822 823
  /* 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
824
      int bid = buff[k].ind;
825 826 827
      if (bid != bucket) {
        struct part part = parts[k];
        struct xpart xpart = xparts[k];
828
        struct cell_buff temp_buff = buff[k];
829 830
        while (bid != bucket) {
          int j = bucket_offset[bid] + bucket_count[bid]++;
Matthieu Schaller's avatar
Matthieu Schaller committed
831
          while (buff[j].ind == bid) {
832 833 834
            j++;
            bucket_count[bid]++;
          }
Pedro Gonnet's avatar
Pedro Gonnet committed
835 836
          memswap(&parts[j], &part, sizeof(struct part));
          memswap(&xparts[j], &xpart, sizeof(struct xpart));
837
          memswap(&buff[j], &temp_buff, sizeof(struct cell_buff));
838 839
          if (parts[j].gpart)
            parts[j].gpart->id_or_neg_offset = -(j + parts_offset);
840
          bid = temp_buff.ind;
841 842 843
        }
        parts[k] = part;
        xparts[k] = xpart;
844
        buff[k] = temp_buff;
845 846
        if (parts[k].gpart)
          parts[k].gpart->id_or_neg_offset = -(k + parts_offset);
847
      }
848
      bucket_count[bid]++;
849 850 851 852
    }
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
853
  for (int k = 0; k < 8; k++) {
854 855 856
    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]];
857 858
  }

859
#ifdef SWIFT_DEBUG_CHECKS
860 861 862 863 864 865 866 867
  /* 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);
  }

868
  /* Verify that _all_ the parts have been assigned to a cell. */
869 870 871 872 873 874 875 876
  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).");
877 878

  /* Verify a few sub-cells. */
879
  for (int k = 0; k < c->progeny[0]->count; k++)
880 881 882
    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])
883 884
      error("Sorting failed (progeny=0).");
  for (int k = 0; k < c->progeny[1]->count; k++)
885 886 887
    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])
888 889
      error("Sorting failed (progeny=1).");
  for (int k = 0; k < c->progeny[2]->count; k++)
890 891 892
    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])
893
      error("Sorting failed (progeny=2).");
894
  for (int k = 0; k < c->progeny[3]->count; k++)
895 896 897
    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])
898 899
      error("Sorting failed (progeny=3).");
  for (int k = 0; k < c->progeny[4]->count; k++)
900 901 902
    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])
903 904
      error("Sorting failed (progeny=4).");
  for (int k = 0; k < c->progeny[5]->count; k++)
905 906 907
    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])
908 909
      error("Sorting failed (progeny=5).");
  for (int k = 0; k < c->progeny[6]->count; k++)
910 911 912
    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])
913 914
      error("Sorting failed (progeny=6).");
  for (int k = 0; k < c->progeny[7]->count; k++)
915 916 917
    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])
918
      error("Sorting failed (progeny=7).");
919
#endif
920

921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954
  /* 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));
955 956
          if (sparts[j].gpart)
            sparts[j].gpart->id_or_neg_offset = -(j + sparts_offset);
957 958 959 960
          bid = temp_buff.ind;
        }
        sparts[k] = spart;
        sbuff[k] = temp_buff;
961 962
        if (sparts[k].gpart)
          sparts[k].gpart->id_or_neg_offset = -(k + sparts_offset);
963 964 965 966 967 968 969 970 971 972 973 974
      }
      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]];
  }

  /* Finally, do the same song and dance for the gparts. */
975 976 977 978
  for (int k = 0; k < 8; k++) bucket_count[k] = 0;

  /* Fill the buffer with the indices. */
  for (int k = 0; k < gcount; k++) {
979 980
    const int bid = (gbuff[k].x[0] > pivot[0]) * 4 +
                    (gbuff[k].x[1] > pivot[1]) * 2 + (gbuff[k].x[2] > pivot[2]);
981
    bucket_count[bid]++;
982
    gbuff[k].ind = bid;
983
  }
984 985 986 987 988 989

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

992 993 994 995
  /* 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++) {
996
      int bid = gbuff[k].ind;
997 998
      if (bid != bucket) {
        struct gpart gpart = gparts[k];
999
        struct cell_buff temp_buff = gbuff[k];
1000 1001
        while (bid != bucket) {
          int j = bucket_offset[bid] + bucket_count[bid]++;
1002
          while (gbuff[j].ind == bid) {
1003 1004 1005
            j++;
            bucket_count[bid]++;
          }
Pedro Gonnet's avatar
Pedro Gonnet committed
1006
          memswap(&gparts[j], &gpart, sizeof(struct gpart));
1007
          memswap(&gbuff[j], &temp_buff, sizeof(struct cell_buff));
1008 1009 1010 1011 1012 1013 1014
          if (gparts[j].type == swift_type_gas) {
            parts[-gparts[j].id_or_neg_offset - parts_offset].gpart =
                &gparts[j];
          } else if (gparts[j].type == swift_type_star) {
            sparts[-gparts[j].id_or_neg_offset - sparts_offset].gpart =
                &gparts[j];
          }
1015
          bid = temp_buff.ind;
1016 1017
        }
        gparts[k] = gpart;
1018
        gbuff[k] = temp_buff;
1019 1020 1021 1022 1023 1024
        if (gparts[k].type == swift_type_gas) {
          parts[-gparts[k].id_or_neg_offset - parts_offset].gpart = &gparts[k];
        } else if (gparts[k].type == swift_type_star) {
          sparts[-gparts[k].id_or_neg_offset - sparts_offset].gpart =
              &gparts[k];
        }
1025
      }
1026
      bucket_count[bid]++;
1027 1028 1029 1030
    }
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
1031
  for (int k = 0; k < 8; k++) {
1032 1033
    c->progeny[k]->gcount = bucket_count[k];
    c->progeny[k]->gparts = &c->gparts[bucket_offset[k]];
1034 1035
  }
}
1036

1037 1038 1039 1040
/**
 * @brief Sanitizes the smoothing length values of cells by setting large
 * outliers to more sensible values.
 *
1041 1042
 * 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.
1043 1044
 *
 * @param c The cell.
1045
 * @param treated Has the cell already been sanitized at this level ?
1046
 */
1047
void cell_sanitize(struct cell *c, int treated) {
1048 1049 1050

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

1053 1054
  /* Treat cells will <1000 particles */
  if (count < 1000 && !treated) {
1055

1056 1057
    /* Get an upper bound on h */
    const float upper_h_max = c->dmin / (1.2f * kernel_gamma);
1058

1059 1060 1061 1062 1063 1064
    /* 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;
    }
  }
1065

1066 1067
  /* Recurse and gather the new h_max values */
  if (c->split) {
1068

1069 1070
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
1071

1072 1073
        /* Recurse */
        cell_sanitize(c->progeny[k], (count < 1000));
1074

1075 1076 1077 1078
        /* And collect */
        h_max = max(h_max, c->progeny[k]->h_max);
      }
    }
1079 1080
  } else {

1081 1082
    /* Get the new value of h_max */
    for (int i = 0; i < count; ++i) h_max = max(h_max, parts[i].h);
1083
  }
1084 1085 1086

  /* Record the change */
  c->h_max = h_max;
1087 1088
}

Matthieu Schaller's avatar
Matthieu Schaller committed
1089 1090 1091 1092 1093 1094
/**
 * @brief Cleans the links in a given cell.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
1095
void cell_clean_links(struct cell *c, void *data) {