cell.c 17.8 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 "atomic.h"
51
#include "error.h"
52
#include "gravity.h"
53
#include "hydro.h"
54
#include "hydro_properties.h"
55
56
#include "space.h"
#include "timers.h"
57

58
59
60
/* Global variables. */
int cell_next_tag = 0;

61
62
63
64
65
66
/**
 * @brief Get the size of the cell subtree.
 *
 * @param c The #cell.
 */

67
int cell_getsize(struct cell *c) {
68

Pedro Gonnet's avatar
Pedro Gonnet committed
69
70
  /* Number of cells in this subtree. */
  int count = 1;
71

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

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

/**
82
83
84
85
86
87
88
89
90
 * @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.
 */

91
92
93
94
int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {

  /* Unpack the current pcell. */
  c->h_max = pc->h_max;
95
96
  c->ti_end_min = pc->ti_end_min;
  c->ti_end_max = pc->ti_end_max;
97
  c->count = pc->count;
98
  c->gcount = pc->gcount;
99
  c->tag = pc->tag;
Matthieu Schaller's avatar
Matthieu Schaller committed
100

101
102
  /* Number of new cells created. */
  int count = 1;
103
104

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
105
  for (int k = 0; k < 8; k++)
106
    if (pc->progeny[k] >= 0) {
Pedro Gonnet's avatar
Pedro Gonnet committed
107
      struct cell *temp = space_getcell(s);
108
      temp->count = 0;
109
      temp->gcount = 0;
110
111
112
113
114
115
116
117
118
119
120
121
      temp->loc[0] = c->loc[0];
      temp->loc[1] = c->loc[1];
      temp->loc[2] = c->loc[2];
      temp->h[0] = c->h[0] / 2;
      temp->h[1] = c->h[1] / 2;
      temp->h[2] = c->h[2] / 2;
      temp->dmin = c->dmin / 2;
      if (k & 4) temp->loc[0] += temp->h[0];
      if (k & 2) temp->loc[1] += temp->h[1];
      if (k & 1) temp->loc[2] += temp->h[2];
      temp->depth = c->depth + 1;
      temp->split = 0;
122
      temp->dx_max = 0.f;
123
124
125
126
127
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
      count += cell_unpack(&pc[pc->progeny[k]], temp, s);
128
129
    }

130
  /* Return the total number of unpacked cells. */
131
  c->pcell_size = count;
132
133
  return count;
}
134

135
/**
136
 * @brief Link the cells recursively to the given #part array.
137
138
139
140
141
142
143
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */

144
int cell_link_parts(struct cell *c, struct part *parts) {
145

146
147
148
  c->parts = parts;

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
149
150
151
152
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
153
        offset += cell_link_parts(c->progeny[k], &parts[offset]);
Pedro Gonnet's avatar
Pedro Gonnet committed
154
155
    }
  }
156

157
  /* Return the total number of linked particles. */
158
159
  return c->count;
}
160

161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
/**
 * @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;
}

187
188
189
190
191
192
193
194
195
196
/**
 * @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.
 */

197
198
199
200
int cell_pack(struct cell *c, struct pcell *pc) {

  /* Start by packing the data of the current cell. */
  pc->h_max = c->h_max;
201
202
  pc->ti_end_min = c->ti_end_min;
  pc->ti_end_max = c->ti_end_max;
203
  pc->count = c->count;
204
  pc->gcount = c->gcount;
205
206
207
  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
208
209
  int count = 1;
  for (int k = 0; k < 8; k++)
210
211
212
213
214
215
216
    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. */
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
  c->pcell_size = count;
  return count;
}

int cell_pack_ti_ends(struct cell *c, int *ti_ends) {

  /* Pack this cell's data. */
  ti_ends[0] = c->ti_end_min;
  
  /* 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;
}

int cell_unpack_ti_ends(struct cell *c, int *ti_ends) {

  /* Unpack this cell's data. */
  c->ti_end_min = ti_ends[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_ti_ends(c->progeny[k], &ti_ends[count]);
    }

  /* Return the number of packed values. */
250
251
  return count;
}
252

253
254
255
256
257
258
/**
 * @brief Lock a cell and hold its parents.
 *
 * @param c The #cell.
 */

259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
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
281
  struct cell *finger;
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
  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
304
305
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
306
      atomic_dec(&finger2->hold);
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338

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

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

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
339
  struct cell *finger;
340
341
342
343
344
345
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

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

    /* Increment the hold. */
346
    atomic_inc(&finger->ghold);
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361

    /* 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
362
363
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
364
      atomic_dec(&finger2->ghold);
365
366
367
368
369
370
371
372
373

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

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

375
/**
376
 * @brief Unlock a cell's parents.
377
378
379
 *
 * @param c The #cell.
 */
380
381
382
383
384
385
386
387
388

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
389
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
390
    atomic_dec(&finger->hold);
391
392
393
394
395
396
397
398
399
400
401
402

  TIMER_TOC(timer_locktree);
}

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
403
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
404
    atomic_dec(&finger->ghold);
405
406
407
408

  TIMER_TOC(timer_locktree);
}

409
410
411
412
413
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
 */
414
415
416

void cell_split(struct cell *c) {

Pedro Gonnet's avatar
Pedro Gonnet committed
417
418
  int i, j;
  const int count = c->count, gcount = c->gcount;
419
420
421
  struct part *parts = c->parts;
  struct xpart *xparts = c->xparts;
  struct gpart *gparts = c->gparts;
422
423
424
425
  int left[8], right[8];
  double pivot[3];

  /* Init the pivots. */
Pedro Gonnet's avatar
Pedro Gonnet committed
426
  for (int k = 0; k < 3; k++) pivot[k] = c->loc[k] + c->h[k] / 2;
427
428
429
430
431
432
433
434

  /* Split along the x-axis. */
  i = 0;
  j = count - 1;
  while (i <= j) {
    while (i <= count - 1 && parts[i].x[0] <= pivot[0]) i += 1;
    while (j >= 0 && parts[j].x[0] > pivot[0]) j -= 1;
    if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
435
      struct part temp = parts[i];
436
437
      parts[i] = parts[j];
      parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
438
      struct xpart xtemp = xparts[i];
439
440
441
442
      xparts[i] = xparts[j];
      xparts[j] = xtemp;
    }
  }
443
444
445
446
447
448
449
450

#ifdef SWIFT_DEBUG_CHECKS
  for (int k = 0; k <= j; k++)
    if (parts[k].x[0] > pivot[0]) error("cell_split: sorting failed.");
  for (int k = i; k < count; k++)
    if (parts[k].x[0] < pivot[0]) error("cell_split: sorting failed.");
#endif

451
452
453
454
455
456
  left[1] = i;
  right[1] = count - 1;
  left[0] = 0;
  right[0] = j;

  /* Split along the y axis, twice. */
Pedro Gonnet's avatar
Pedro Gonnet committed
457
  for (int k = 1; k >= 0; k--) {
458
459
460
461
462
463
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && parts[i].x[1] <= pivot[1]) i += 1;
      while (j >= left[k] && parts[j].x[1] > pivot[1]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
464
        struct part temp = parts[i];
465
466
        parts[i] = parts[j];
        parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
467
        struct xpart xtemp = xparts[i];
468
469
470
471
        xparts[i] = xparts[j];
        xparts[j] = xtemp;
      }
    }
472
473
474
475
476
477
478
479
480
481
482

#ifdef SWIFT_DEBUG_CHECKS
    for (int kk = left[k]; kk <= j; kk++)
      if (parts[kk].x[1] > pivot[1]) {
        message("ival=[%i,%i], i=%i, j=%i.", left[k], right[k], i, j);
        error("sorting failed (left).");
      }
    for (int kk = i; kk <= right[k]; kk++)
      if (parts[kk].x[1] < pivot[1]) error("sorting failed (right).");
#endif

483
484
485
486
487
488
489
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Split along the z axis, four times. */
Pedro Gonnet's avatar
Pedro Gonnet committed
490
  for (int k = 3; k >= 0; k--) {
491
492
493
494
495
496
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && parts[i].x[2] <= pivot[2]) i += 1;
      while (j >= left[k] && parts[j].x[2] > pivot[2]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
497
        struct part temp = parts[i];
498
499
        parts[i] = parts[j];
        parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
500
        struct xpart xtemp = xparts[i];
501
502
503
504
        xparts[i] = xparts[j];
        xparts[j] = xtemp;
      }
    }
505
506
507
508
509
510
511
512
513
514
515
516
517
518

#ifdef SWIFT_DEBUG_CHECKS
    for (int kk = left[k]; kk <= j; kk++)
      if (parts[kk].x[2] > pivot[2]) {
        message("ival=[%i,%i], i=%i, j=%i.", left[k], right[k], i, j);
        error("sorting failed (left).");
      }
    for (int kk = i; kk <= right[k]; kk++)
      if (parts[kk].x[2] < pivot[2]) {
        message("ival=[%i,%i], i=%i, j=%i.", left[k], right[k], i, j);
        error("sorting failed (right).");
      }
#endif

519
520
521
522
523
524
525
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
526
  for (int k = 0; k < 8; k++) {
527
528
529
530
531
532
    c->progeny[k]->count = right[k] - left[k] + 1;
    c->progeny[k]->parts = &c->parts[left[k]];
    c->progeny[k]->xparts = &c->xparts[left[k]];
  }

  /* Re-link the gparts. */
Pedro Gonnet's avatar
Pedro Gonnet committed
533
  for (int k = 0; k < count; k++)
534
535
    if (parts[k].gpart != NULL) parts[k].gpart->part = &parts[k];

536
#ifdef SWIFT_DEBUG_CHECKS
537
  /* Verify that _all_ the parts have been assigned to a cell. */
538
539
540
541
542
543
544
545
  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).");
546
547

  /* Verify a few sub-cells. */
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
  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).");
#endif
564
565
566
567
568
569
570
571
572
573

  /* Now do the same song and dance for the gparts. */

  /* Split along the x-axis. */
  i = 0;
  j = gcount - 1;
  while (i <= j) {
    while (i <= gcount - 1 && gparts[i].x[0] <= pivot[0]) i += 1;
    while (j >= 0 && gparts[j].x[0] > pivot[0]) j -= 1;
    if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
574
      struct gpart gtemp = gparts[i];
575
576
      gparts[i] = gparts[j];
      gparts[j] = gtemp;
577
    }
578
579
580
581
582
583
584
  }
  left[1] = i;
  right[1] = gcount - 1;
  left[0] = 0;
  right[0] = j;

  /* Split along the y axis, twice. */
Pedro Gonnet's avatar
Pedro Gonnet committed
585
  for (int k = 1; k >= 0; k--) {
586
587
588
589
590
591
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && gparts[i].x[1] <= pivot[1]) i += 1;
      while (j >= left[k] && gparts[j].x[1] > pivot[1]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
592
        struct gpart gtemp = gparts[i];
593
594
595
596
597
598
599
600
601
602
603
        gparts[i] = gparts[j];
        gparts[j] = gtemp;
      }
    }
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Split along the z axis, four times. */
Pedro Gonnet's avatar
Pedro Gonnet committed
604
  for (int k = 3; k >= 0; k--) {
605
606
607
608
609
610
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && gparts[i].x[2] <= pivot[2]) i += 1;
      while (j >= left[k] && gparts[j].x[2] > pivot[2]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
611
        struct gpart gtemp = gparts[i];
612
613
614
615
616
617
618
619
620
621
622
        gparts[i] = gparts[j];
        gparts[j] = gtemp;
      }
    }
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
623
  for (int k = 0; k < 8; k++) {
624
625
626
627
628
    c->progeny[k]->gcount = right[k] - left[k] + 1;
    c->progeny[k]->gparts = &c->gparts[left[k]];
  }

  /* Re-link the parts. */
Pedro Gonnet's avatar
Pedro Gonnet committed
629
  for (int k = 0; k < gcount; k++)
630
631
    if (gparts[k].id > 0) gparts[k].part->gpart = &gparts[k];
}
632
633
634
635
636
637
638
639
640
641
642

/**
 * @brief Initialises all particles to a valid state even if the ICs were stupid
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_init_parts(struct cell *c, void *data) {

  struct part *p = c->parts;
  struct xpart *xp = c->xparts;
643
  const int count = c->count;
644

645
  for (int i = 0; i < count; ++i) {
646
647
    p[i].ti_begin = 0;
    p[i].ti_end = 0;
648
649
650
    xp[i].v_full[0] = p[i].v[0];
    xp[i].v_full[1] = p[i].v[1];
    xp[i].v_full[2] = p[i].v[2];
651
    hydro_first_init_part(&p[i], &xp[i]);
652
653
    hydro_init_part(&p[i]);
    hydro_reset_acceleration(&p[i]);
654
  }
655
656
  c->ti_end_min = 0;
  c->ti_end_max = 0;
657
658
}

659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
/**
 * @brief Initialises all g-particles to a valid state even if the ICs were
 *stupid
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_init_gparts(struct cell *c, void *data) {

  struct gpart *gp = c->gparts;
  const int gcount = c->gcount;

  for (int i = 0; i < gcount; ++i) {
    gp[i].ti_begin = 0;
    gp[i].ti_end = 0;
674
    gravity_first_init_gpart(&gp[i]);
675
676
677
678
679
  }
  c->ti_end_min = 0;
  c->ti_end_max = 0;
}

680
/**
681
682
 * @brief Converts hydro quantities to a valid state after the initial density
 *calculation
683
684
685
686
687
688
689
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_convert_hydro(struct cell *c, void *data) {

  struct part *p = c->parts;
690
691

  for (int i = 0; i < c->count; ++i) {
692
693
694
695
    hydro_convert_quantities(&p[i]);
  }
}

Matthieu Schaller's avatar
Matthieu Schaller committed
696
697
698
699
700
701
/**
 * @brief Cleans the links in a given cell.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
702
void cell_clean_links(struct cell *c, void *data) {
Matthieu Schaller's avatar
Matthieu Schaller committed
703
704
  c->density = NULL;
  c->nr_density = 0;
705

Matthieu Schaller's avatar
Matthieu Schaller committed
706
707
708
  c->force = NULL;
  c->nr_force = 0;
}