cell.c 17.1 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
55
#include "space.h"
#include "timers.h"
56

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

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

66
int cell_getsize(struct cell *c) {
67

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

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

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

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

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

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

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

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

129
130
131
  /* Return the total number of unpacked cells. */
  return count;
}
132

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

142
int cell_link_parts(struct cell *c, struct part *parts) {
143

144
145
146
  c->parts = parts;

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

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

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

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

195
196
197
198
int cell_pack(struct cell *c, struct pcell *pc) {

  /* Start by packing the data of the current cell. */
  pc->h_max = c->h_max;
199
200
  pc->ti_end_min = c->ti_end_min;
  pc->ti_end_max = c->ti_end_max;
201
  pc->count = c->count;
202
  pc->gcount = c->gcount;
203
204
205
  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
206
207
  int count = 1;
  for (int k = 0; k < 8; k++)
208
209
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. */
  return count;
}
217

218
219
220
221
222
223
/**
 * @brief Lock a cell and hold its parents.
 *
 * @param c The #cell.
 */

224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
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
246
  struct cell *finger;
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
  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
269
270
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
      __sync_fetch_and_sub(&finger2->hold, 1);

    /* 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
304
  struct cell *finger;
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

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

    /* Increment the hold. */
    __sync_fetch_and_add(&finger->ghold, 1);

    /* 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
327
328
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
329
330
331
332
333
334
335
336
337
338
      __sync_fetch_and_sub(&finger2->ghold, 1);

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

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

340
/**
341
 * @brief Unlock a cell's parents.
342
343
344
 *
 * @param c The #cell.
 */
345
346
347
348
349
350
351
352
353

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
354
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
355
356
357
358
359
360
361
362
363
364
365
366
367
    __sync_fetch_and_sub(&finger->hold, 1);

  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
368
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
369
370
371
372
373
    __sync_fetch_and_sub(&finger->ghold, 1);

  TIMER_TOC(timer_locktree);
}

374
375
376
377
378
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
 */
379
380
381

void cell_split(struct cell *c) {

Pedro Gonnet's avatar
Pedro Gonnet committed
382
383
  int i, j;
  const int count = c->count, gcount = c->gcount;
384
385
386
  struct part *parts = c->parts;
  struct xpart *xparts = c->xparts;
  struct gpart *gparts = c->gparts;
387
388
389
390
  int left[8], right[8];
  double pivot[3];

  /* Init the pivots. */
Pedro Gonnet's avatar
Pedro Gonnet committed
391
  for (int k = 0; k < 3; k++) pivot[k] = c->loc[k] + c->h[k] / 2;
392
393
394
395
396
397
398
399

  /* 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
400
      struct part temp = parts[i];
401
402
      parts[i] = parts[j];
      parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
403
      struct xpart xtemp = xparts[i];
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
      xparts[i] = xparts[j];
      xparts[j] = xtemp;
    }
  }
  /* for ( k = 0 ; k <= j ; k++ )
      if ( parts[k].x[0] > pivot[0] )
          error( "cell_split: sorting failed." );
  for ( k = i ; k < count ; k++ )
      if ( parts[k].x[0] < pivot[0] )
          error( "cell_split: sorting failed." ); */
  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
420
  for (int k = 1; k >= 0; k--) {
421
422
423
424
425
426
    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
427
        struct part temp = parts[i];
428
429
        parts[i] = parts[j];
        parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
430
        struct xpart xtemp = xparts[i];
431
432
433
434
435
436
437
438
        xparts[i] = xparts[j];
        xparts[j] = xtemp;
      }
    }
    /* 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)." );
439
            }
440
441
442
443
444
445
446
447
448
449
    for ( int kk = i ; kk <= right[k] ; kk++ )
        if ( parts[kk].x[1] < pivot[1] )
            error( "sorting failed (right)." ); */
    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
450
  for (int k = 3; k >= 0; k--) {
451
452
453
454
455
456
    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
457
        struct part temp = parts[i];
458
459
        parts[i] = parts[j];
        parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
460
        struct xpart xtemp = xparts[i];
461
462
463
464
465
466
467
468
        xparts[i] = xparts[j];
        xparts[j] = xtemp;
      }
    }
    /* 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)." );
469
            }
470
471
472
473
474
475
476
477
478
479
480
481
    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)." );
            } */
    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
482
  for (int k = 0; k < 8; k++) {
483
484
485
486
487
488
    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
489
  for (int k = 0; k < count; k++)
490
491
492
493
494
495
496
497
498
499
500
501
502
    if (parts[k].gpart != NULL) parts[k].gpart->part = &parts[k];

  /* Verify that _all_ the parts have been assigned to a cell. */
  /* for ( 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)." ); */

  /* Verify a few sub-cells. */
Pedro Gonnet's avatar
Pedro Gonnet committed
503
  /* for (int k = 0 ; k < c->progeny[0]->count ; k++ )
504
505
506
507
      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)." );
Pedro Gonnet's avatar
Pedro Gonnet committed
508
  for (int k = 0 ; k < c->progeny[1]->count ; k++ )
509
510
511
512
      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)." );
Pedro Gonnet's avatar
Pedro Gonnet committed
513
  for (int k = 0 ; k < c->progeny[2]->count ; k++ )
514
515
516
517
518
519
520
521
522
523
524
525
526
527
      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)." ); */

  /* 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
528
      struct gpart gtemp = gparts[i];
529
530
      gparts[i] = gparts[j];
      gparts[j] = gtemp;
531
    }
532
533
534
535
536
537
538
  }
  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
539
  for (int k = 1; k >= 0; k--) {
540
541
542
543
544
545
    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
546
        struct gpart gtemp = gparts[i];
547
548
549
550
551
552
553
554
555
556
557
        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
558
  for (int k = 3; k >= 0; k--) {
559
560
561
562
563
564
    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
565
        struct gpart gtemp = gparts[i];
566
567
568
569
570
571
572
573
574
575
576
        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
577
  for (int k = 0; k < 8; k++) {
578
579
580
581
582
    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
583
  for (int k = 0; k < gcount; k++)
584
585
    if (gparts[k].id > 0) gparts[k].part->gpart = &gparts[k];
}
586
587
588
589
590
591
592
593
594
595
596

/**
 * @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;
597
  const int count = c->count;
598

599
  for (int i = 0; i < count; ++i) {
600
601
    p[i].ti_begin = 0;
    p[i].ti_end = 0;
602
603
604
    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];
605
    hydro_first_init_part(&p[i], &xp[i]);
606
607
    hydro_init_part(&p[i]);
    hydro_reset_acceleration(&p[i]);
608
  }
609
610
  c->ti_end_min = 0;
  c->ti_end_max = 0;
611
612
}

613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
/**
 * @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;
628
    gravity_first_init_gpart(&gp[i]);
629
630
631
632
633
  }
  c->ti_end_min = 0;
  c->ti_end_max = 0;
}

634
/**
635
636
 * @brief Converts hydro quantities to a valid state after the initial density
 *calculation
637
638
639
640
641
642
643
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_convert_hydro(struct cell *c, void *data) {

  struct part *p = c->parts;
644
645

  for (int i = 0; i < c->count; ++i) {
646
647
648
649
    hydro_convert_quantities(&p[i]);
  }
}

Matthieu Schaller's avatar
Matthieu Schaller committed
650
651
652
653
654
655
/**
 * @brief Cleans the links in a given cell.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
656
void cell_clean_links(struct cell *c, void *data) {
Matthieu Schaller's avatar
Matthieu Schaller committed
657
658
  c->density = NULL;
  c->nr_density = 0;
659

Matthieu Schaller's avatar
Matthieu Schaller committed
660
661
662
  c->force = NULL;
  c->nr_force = 0;
}