runner.c 53.9 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
 ******************************************************************************/
Pedro Gonnet's avatar
Pedro Gonnet committed
23

Pedro Gonnet's avatar
Pedro Gonnet committed
24
25
/* Config parameters. */
#include "../config.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
26
27
28
29

/* Some standard headers. */
#include <float.h>
#include <limits.h>
30
#include <stdlib.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
31

32
33
/* MPI headers. */
#ifdef WITH_MPI
34
#include <mpi.h>
35
36
#endif

37
38
39
/* This object's header. */
#include "runner.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
40
/* Local headers. */
41
#include "active.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
42
#include "approx_math.h"
43
#include "atomic.h"
44
#include "cell.h"
45
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
46
#include "cooling.h"
47
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
48
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
49
#include "engine.h"
50
#include "error.h"
51
52
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
53
#include "hydro_properties.h"
54
#include "kick.h"
55
#include "minmax.h"
James Willis's avatar
James Willis committed
56
#include "runner_doiact_vec.h"
57
#include "scheduler.h"
58
#include "sourceterms.h"
59
#include "space.h"
60
#include "stars.h"
61
62
#include "task.h"
#include "timers.h"
63
#include "timestep.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
64

65
66
67
68
69
70
71
72
73
74
75
76
/**
 * @brief  Entry in a list of sorted indices.
 */
struct entry {

  /*! Distance on the axis */
  float d;

  /*! Particle index */
  int i;
};

77
/* Orientation of the cell pairs */
Matthieu Schaller's avatar
Matthieu Schaller committed
78
79
const double runner_shift[13][3] = {
    {5.773502691896258e-01, 5.773502691896258e-01, 5.773502691896258e-01},
80
    {7.071067811865475e-01, 7.071067811865475e-01, 0.0},
Matthieu Schaller's avatar
Matthieu Schaller committed
81
82
83
84
85
86
87
88
89
90
    {5.773502691896258e-01, 5.773502691896258e-01, -5.773502691896258e-01},
    {7.071067811865475e-01, 0.0, 7.071067811865475e-01},
    {1.0, 0.0, 0.0},
    {7.071067811865475e-01, 0.0, -7.071067811865475e-01},
    {5.773502691896258e-01, -5.773502691896258e-01, 5.773502691896258e-01},
    {7.071067811865475e-01, -7.071067811865475e-01, 0.0},
    {5.773502691896258e-01, -5.773502691896258e-01, -5.773502691896258e-01},
    {0.0, 7.071067811865475e-01, 7.071067811865475e-01},
    {0.0, 1.0, 0.0},
    {0.0, 7.071067811865475e-01, -7.071067811865475e-01},
91
92
    {0.0, 0.0, 1.0},
};
93
94

/* Does the axis need flipping ? */
95
96
const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
                              0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
Pedro Gonnet's avatar
Pedro Gonnet committed
97

98
/* Import the density loop functions. */
99
100
101
#define FUNCTION density
#include "runner_doiact.h"

102
/* Import the gradient loop functions (if required). */
103
104
105
106
107
108
#ifdef EXTRA_HYDRO_LOOP
#undef FUNCTION
#define FUNCTION gradient
#include "runner_doiact.h"
#endif

109
/* Import the force loop functions. */
110
111
112
113
#undef FUNCTION
#define FUNCTION force
#include "runner_doiact.h"

114
/* Import the gravity loop functions. */
115
#include "runner_doiact_fft.h"
116
#include "runner_doiact_grav.h"
117

Tom Theuns's avatar
Tom Theuns committed
118
/**
Tom Theuns's avatar
Tom Theuns committed
119
 * @brief Perform source terms
Tom Theuns's avatar
Tom Theuns committed
120
121
122
123
124
125
126
 *
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
 */
void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
  const int count = c->count;
127
  const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
Tom Theuns's avatar
Tom Theuns committed
128
  const double cell_width[3] = {c->width[0], c->width[1], c->width[2]};
Tom Theuns's avatar
Tom Theuns committed
129
  struct sourceterms *sourceterms = r->e->sourceterms;
130
  const int dimen = 3;
Tom Theuns's avatar
Tom Theuns committed
131
132
133
134
135
136
137

  TIMER_TIC;

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_sourceterms(r, c->progeny[k], 0);
138
  } else {
Tom Theuns's avatar
Tom Theuns committed
139

140
    if (count > 0) {
Tom Theuns's avatar
Tom Theuns committed
141

142
143
144
145
146
147
      /* do sourceterms in this cell? */
      const int incell =
          sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
      if (incell == 1) {
        sourceterms_apply(r, sourceterms, c);
      }
Tom Theuns's avatar
Tom Theuns committed
148
149
    }
  }
Tom Theuns's avatar
Tom Theuns committed
150
151
152
153

  if (timer) TIMER_TOC(timer_dosource);
}

Tom Theuns's avatar
Tom Theuns committed
154
155
156
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
157
158
159
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
160
 */
161
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
162

Matthieu Schaller's avatar
Matthieu Schaller committed
163
164
  struct gpart *restrict gparts = c->gparts;
  const int gcount = c->gcount;
165
166
167
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
168
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
169

170
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
171

172
  /* Anything to do here? */
173
  if (!cell_is_active(c, e)) return;
174

Tom Theuns's avatar
Tom Theuns committed
175
176
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
177
    for (int k = 0; k < 8; k++)
178
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
179
  } else {
180

181
182
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
183

184
185
      /* Get a direct pointer on the part. */
      struct gpart *restrict gp = &gparts[i];
Matthieu Schaller's avatar
Matthieu Schaller committed
186

187
      /* Is this part within the time step? */
188
      if (gpart_is_active(gp, e)) {
189
190
        external_gravity_acceleration(time, potential, constants, gp);
      }
191
    }
192
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
193

194
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
195
196
}

Stefan Arridge's avatar
Stefan Arridge committed
197
/**
198
199
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
200
201
202
203
204
205
206
207
 *
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
 */
void runner_do_cooling(struct runner *r, struct cell *c, int timer) {

  struct part *restrict parts = c->parts;
208
  struct xpart *restrict xparts = c->xparts;
Stefan Arridge's avatar
Stefan Arridge committed
209
  const int count = c->count;
210
211
212
  const struct engine *e = r->e;
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
213
  const struct unit_system *us = e->internal_units;
214
  const double timeBase = e->timeBase;
Stefan Arridge's avatar
Stefan Arridge committed
215
216
217

  TIMER_TIC;

218
219
220
  /* Anything to do here? */
  if (!cell_is_active(c, e)) return;

Stefan Arridge's avatar
Stefan Arridge committed
221
222
223
224
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
225
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
226

227
228
    /* Loop over the parts in this cell. */
    for (int i = 0; i < count; i++) {
Stefan Arridge's avatar
Stefan Arridge committed
229

230
231
232
      /* Get a direct pointer on the part. */
      struct part *restrict p = &parts[i];
      struct xpart *restrict xp = &xparts[i];
Stefan Arridge's avatar
Stefan Arridge committed
233

234
      if (part_is_active(p, e)) {
235

236
237
        /* Let's cool ! */
        const double dt = get_timestep(p->time_bin, timeBase);
238
239
        cooling_cool_part(constants, us, cooling_func, p, xp, dt);
      }
Stefan Arridge's avatar
Stefan Arridge committed
240
241
242
243
244
245
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
246
247
248
249
250
251
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
252
void runner_do_sort_ascending(struct entry *sort, int N) {
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
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306

  struct {
    short int lo, hi;
  } qstack[10];
  int qpos, i, j, lo, hi, imin;
  struct entry temp;
  float pivot;

  /* Sort parts in cell_i in decreasing order with quicksort */
  qstack[0].lo = 0;
  qstack[0].hi = N - 1;
  qpos = 0;
  while (qpos >= 0) {
    lo = qstack[qpos].lo;
    hi = qstack[qpos].hi;
    qpos -= 1;
    if (hi - lo < 15) {
      for (i = lo; i < hi; i++) {
        imin = i;
        for (j = i + 1; j <= hi; j++)
          if (sort[j].d < sort[imin].d) imin = j;
        if (imin != i) {
          temp = sort[imin];
          sort[imin] = sort[i];
          sort[i] = temp;
        }
      }
    } else {
      pivot = sort[(lo + hi) / 2].d;
      i = lo;
      j = hi;
      while (i <= j) {
        while (sort[i].d < pivot) i++;
        while (sort[j].d > pivot) j--;
        if (i <= j) {
          if (i < j) {
            temp = sort[i];
            sort[i] = sort[j];
            sort[j] = temp;
          }
          i += 1;
          j -= 1;
        }
      }
      if (j > (lo + hi) / 2) {
        if (lo < j) {
          qpos += 1;
          qstack[qpos].lo = lo;
          qstack[qpos].hi = j;
        }
        if (i < hi) {
          qpos += 1;
          qstack[qpos].lo = i;
          qstack[qpos].hi = hi;
Pedro Gonnet's avatar
Pedro Gonnet committed
307
        }
308
309
310
311
312
313
314
315
316
317
318
319
      } else {
        if (i < hi) {
          qpos += 1;
          qstack[qpos].lo = i;
          qstack[qpos].hi = hi;
        }
        if (lo < j) {
          qpos += 1;
          qstack[qpos].lo = lo;
          qstack[qpos].hi = j;
        }
      }
Pedro Gonnet's avatar
Pedro Gonnet committed
320
    }
321
322
323
  }
}

Pedro Gonnet's avatar
Pedro Gonnet committed
324
325
326
327
328
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
329
 * @param flags Cell flag.
330
331
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
332
 */
333
void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
334
335
336
337
338
339
340

  struct entry *finger;
  struct entry *fingers[8];
  struct part *parts = c->parts;
  struct entry *sort;
  int j, k, count = c->count;
  int i, ind, off[8], inds[8], temp_i, missing;
Matthieu Schaller's avatar
Matthieu Schaller committed
341
342
  float buff[8];
  double px[3];
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366

  TIMER_TIC

  /* Clean-up the flags, i.e. filter out what's already been sorted. */
  flags &= ~c->sorted;
  if (flags == 0) return;

  /* start by allocating the entry arrays. */
  if (c->sort == NULL || c->sortsize < count) {
    if (c->sort != NULL) free(c->sort);
    c->sortsize = count * 1.1;
    if ((c->sort = (struct entry *)malloc(sizeof(struct entry) *
                                          (c->sortsize + 1) * 13)) == NULL)
      error("Failed to allocate sort memory.");
  }
  sort = c->sort;

  /* Does this cell have any progeny? */
  if (c->split) {

    /* Fill in the gaps within the progeny. */
    for (k = 0; k < 8; k++) {
      if (c->progeny[k] == NULL) continue;
      missing = flags & ~c->progeny[k]->sorted;
367
      if (missing) runner_do_sort(r, c->progeny[k], missing, 0);
368
369
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
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
    }

    /* Loop over the 13 different sort arrays. */
    for (j = 0; j < 13; j++) {

      /* Has this sort array been flagged? */
      if (!(flags & (1 << j))) continue;

      /* Init the particle index offsets. */
      for (off[0] = 0, k = 1; k < 8; k++)
        if (c->progeny[k - 1] != NULL)
          off[k] = off[k - 1] + c->progeny[k - 1]->count;
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
      for (k = 0; k < 8; k++) {
        inds[k] = k;
        if (c->progeny[k] != NULL && c->progeny[k]->count > 0) {
          fingers[k] = &c->progeny[k]->sort[j * (c->progeny[k]->count + 1)];
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
      for (i = 0; i < 7; i++)
        for (k = i + 1; k < 8; k++)
          if (buff[inds[k]] < buff[inds[i]]) {
            temp_i = inds[i];
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
      finger = &sort[j * (count + 1)];
      for (ind = 0; ind < count; ind++) {

        /* Copy the minimum into the new sort array. */
        finger[ind].d = buff[inds[0]];
        finger[ind].i = fingers[inds[0]]->i + off[inds[0]];

        /* Update the buffer. */
        fingers[inds[0]] += 1;
        buff[inds[0]] = fingers[inds[0]]->d;

        /* Find the smallest entry. */
        for (k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
          temp_i = inds[k - 1];
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
420
        }
421

422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
      } /* Merge. */

      /* Add a sentinel. */
      sort[j * (count + 1) + count].d = FLT_MAX;
      sort[j * (count + 1) + count].i = 0;

      /* Mark as sorted. */
      c->sorted |= (1 << j);

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

    /* Fill the sort array. */
    for (k = 0; k < count; k++) {
      px[0] = parts[k].x[0];
      px[1] = parts[k].x[1];
      px[2] = parts[k].x[2];
      for (j = 0; j < 13; j++)
        if (flags & (1 << j)) {
          sort[j * (count + 1) + k].i = k;
Matthieu Schaller's avatar
Matthieu Schaller committed
446
447
448
          sort[j * (count + 1) + k].d = px[0] * runner_shift[j][0] +
                                        px[1] * runner_shift[j][1] +
                                        px[2] * runner_shift[j][2];
449
        }
450
    }
451
452
453
454
455
456

    /* Add the sentinel and sort. */
    for (j = 0; j < 13; j++)
      if (flags & (1 << j)) {
        sort[j * (count + 1) + count].d = FLT_MAX;
        sort[j * (count + 1) + count].i = 0;
457
        runner_do_sort_ascending(&sort[j * (count + 1)], count);
458
459
460
461
        c->sorted |= (1 << j);
      }
  }

462
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
463
  /* Verify the sorting. */
464
465
466
467
468
469
470
471
472
473
  for (j = 0; j < 13; j++) {
    if (!(flags & (1 << j))) continue;
    finger = &sort[j * (count + 1)];
    for (k = 1; k < count; k++) {
      if (finger[k].d < finger[k - 1].d)
        error("Sorting failed, ascending array.");
      if (finger[k].i >= count) error("Sorting failed, indices borked.");
    }
  }
#endif
474
475
476
477

  if (clock) TIMER_TOC(timer_dosort);
}

478
479
480
481
482
/**
 * @brief Initialize the particles before the density calculation
 *
 * @param r The runner thread.
 * @param c The cell.
Matthieu Schaller's avatar
Matthieu Schaller committed
483
 * @param timer 1 if the time is to be recorded.
484
 */
485
void runner_do_init(struct runner *r, struct cell *c, int timer) {
486

Matthieu Schaller's avatar
Matthieu Schaller committed
487
488
  struct part *restrict parts = c->parts;
  struct gpart *restrict gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
489
  const int count = c->count;
490
  const int gcount = c->gcount;
491
  const struct engine *e = r->e;
492
  const struct space *s = e->s;
493
494

  TIMER_TIC;
495

496
  /* Anything to do here? */
497
  if (!cell_is_active(c, e)) return;
498

499
  /* Reset the gravity acceleration tensors */
500
  if (e->policy & engine_policy_self_gravity)
501
    gravity_field_tensors_init(&c->multipole->pot);
502

503
504
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
505
    for (int k = 0; k < 8; k++)
506
      if (c->progeny[k] != NULL) runner_do_init(r, c->progeny[k], 0);
507
508
  } else {

509
510
    /* Loop over the parts in this cell. */
    for (int i = 0; i < count; i++) {
511

512
      /* Get a direct pointer on the part. */
513
      struct part *restrict p = &parts[i];
514

515
      if (part_is_active(p, e)) {
Matthieu Schaller's avatar
Matthieu Schaller committed
516

517
        /* Get ready for a density calculation */
518
        hydro_init_part(p, &s->hs);
519
      }
520
    }
521
522
523
524
525

    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {

      /* Get a direct pointer on the part. */
526
      struct gpart *restrict gp = &gparts[i];
527

528
      if (gpart_is_active(gp, e)) {
529
530

        /* Get ready for a density calculation */
531
        gravity_init_gpart(gp);
532
533
      }
    }
534
  }
535

Peter W. Draper's avatar
Peter W. Draper committed
536
  if (timer) TIMER_TOC(timer_init);
537
538
}

539
/**
540
541
542
543
544
 * @brief Intermediate task after the gradient loop that does final operations
 * on the gradient quantities and optionally slope limits the gradients
 *
 * @param r The runner thread.
 * @param c The cell.
545
 * @param timer Are we timing this ?
546
 */
547
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
548

549
#ifdef EXTRA_HYDRO_LOOP
550

551
552
  struct part *restrict parts = c->parts;
  const int count = c->count;
553
  const struct engine *e = r->e;
554

555
556
  TIMER_TIC;

557
  /* Anything to do here? */
558
  if (!cell_is_active(c, e)) return;
559

560
561
562
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
563
      if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k], 0);
564
565
566
567
568
569
570
571
  } else {

    /* Loop over the parts in this cell. */
    for (int i = 0; i < count; i++) {

      /* Get a direct pointer on the part. */
      struct part *restrict p = &parts[i];

572
      if (part_is_active(p, e)) {
573
574
575
576
577
578

        /* Get ready for a force calculation */
        hydro_end_gradient(p);
      }
    }
  }
579

580
581
  if (timer) TIMER_TOC(timer_do_extra_ghost);

582
583
#else
  error("SWIFT was not compiled with the extra hydro loop activated.");
584
#endif
585
}
586

587
/**
588
589
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
590
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
591
 * @param r The runner thread.
592
 * @param c The cell.
593
 * @param timer Are we timing this ?
594
 */
595
void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
596

Matthieu Schaller's avatar
Matthieu Schaller committed
597
598
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
599
  const struct engine *e = r->e;
600
  const struct space *s = e->s;
601
  const float hydro_h_max = e->hydro_properties->h_max;
602
  const float target_wcount = e->hydro_properties->target_neighbours;
603
  const float max_wcount =
604
      target_wcount + e->hydro_properties->delta_neighbours;
605
  const float min_wcount =
606
607
      target_wcount - e->hydro_properties->delta_neighbours;
  const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
608
  int redo = 0, count = 0;
609

610
611
  TIMER_TIC;

612
  /* Anything to do here? */
613
  if (!cell_is_active(c, e)) return;
614

615
616
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
617
    for (int k = 0; k < 8; k++)
618
619
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k], 0);
  } else {
620

621
    /* Init the list of active particles that have to be updated. */
622
    int *pid = NULL;
623
    if ((pid = malloc(sizeof(int) * c->count)) == NULL)
624
      error("Can't allocate memory for pid.");
625
626
627
628
629
    for (int k = 0; k < c->count; k++)
      if (part_is_active(&parts[k], e)) {
        pid[count] = k;
        ++count;
      }
630

631
632
633
    /* While there are particles that need to be updated... */
    for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
         num_reruns++) {
634

635
636
      /* Reset the redo-count. */
      redo = 0;
637

638
      /* Loop over the remaining active parts in this cell. */
639
      for (int i = 0; i < count; i++) {
640

641
642
643
        /* Get a direct pointer on the part. */
        struct part *restrict p = &parts[pid[i]];
        struct xpart *restrict xp = &xparts[pid[i]];
644

645
#ifdef SWIFT_DEBUG_CHECKS
646
        /* Is this part within the timestep? */
647
648
649
650
651
        if (!part_is_active(p, e)) error("Ghost applied to inactive particle");
#endif

        /* Finish the density calculation */
        hydro_end_density(p);
652

653
654
        /* Did we get the right number of neighbours? */
        if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
655

656
          float h_corr = 0.f;
Matthieu Schaller's avatar
Matthieu Schaller committed
657

658
659
          /* If no derivative, double the smoothing length. */
          if (p->density.wcount_dh == 0.0f) h_corr = p->h;
660

661
662
663
          /* Otherwise, compute the smoothing length update (Newton step). */
          else {
            h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
664

665
666
667
668
            /* Truncate to the range [ -p->h/2 , p->h ]. */
            h_corr = (h_corr < p->h) ? h_corr : p->h;
            h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
          }
669

670
671
          /* Ok, correct then */
          p->h += h_corr;
672

673
674
          /* If below the absolute maximum, try again */
          if (p->h < hydro_h_max) {
675

676
677
678
            /* Flag for another round of fun */
            pid[redo] = pid[i];
            redo += 1;
679

680
            /* Re-initialise everything */
681
            hydro_init_part(p, &s->hs);
682
683
684
685
686
687
688
689

            /* Off we go ! */
            continue;
          } else {

            /* Ok, this particle is a lost cause... */
            p->h = hydro_h_max;
          }
690
        }
691

692
        /* We now have a particle whose smoothing length has converged */
Matthieu Schaller's avatar
Matthieu Schaller committed
693

694
        /* As of here, particle force variables will be set. */
695

696
697
        /* Compute variables required for the force loop */
        hydro_prepare_force(p, xp);
698

699
700
        /* The particle force values are now set.  Do _NOT_
           try to read any particle density variables! */
Matthieu Schaller's avatar
Matthieu Schaller committed
701

702
703
        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
704
705
      }

706
707
      /* We now need to treat the particles whose smoothing length had not
       * converged again */
708

709
710
711
      /* Re-set the counter for the next loop (potentially). */
      count = redo;
      if (count > 0) {
712

713
714
        /* Climb up the cell hierarchy. */
        for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
Matthieu Schaller's avatar
Matthieu Schaller committed
715

716
717
          /* Run through this cell's density interactions. */
          for (struct link *l = finger->density; l != NULL; l = l->next) {
Matthieu Schaller's avatar
Matthieu Schaller committed
718

719
720
721
            /* Self-interaction? */
            if (l->t->type == task_type_self)
              runner_doself_subset_density(r, finger, parts, pid, count);
722

723
724
            /* Otherwise, pair interaction? */
            else if (l->t->type == task_type_pair) {
725

726
727
728
729
730
731
732
              /* Left or right? */
              if (l->t->ci == finger)
                runner_dopair_subset_density(r, finger, parts, pid, count,
                                             l->t->cj);
              else
                runner_dopair_subset_density(r, finger, parts, pid, count,
                                             l->t->ci);
733

734
            }
735

736
737
738
739
            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
              runner_dosub_subset_density(r, finger, parts, pid, count, NULL,
                                          -1, 1);
740

741
742
743
744
745
746
747
748
749
750
751
            /* Otherwise, sub-pair interaction? */
            else if (l->t->type == task_type_sub_pair) {

              /* Left or right? */
              if (l->t->ci == finger)
                runner_dosub_subset_density(r, finger, parts, pid, count,
                                            l->t->cj, -1, 1);
              else
                runner_dosub_subset_density(r, finger, parts, pid, count,
                                            l->t->ci, -1, 1);
            }
752
753
754
          }
        }
      }
755
    }
756

757
758
759
760
761
762
763
#ifdef SWIFT_DEBUG_CHECKS
    if (count) {
      message("Smoothing length failed to converge on %i particles.", count);

      error("Aborting....");
    }
#else
764
765
    if (count)
      message("Smoothing length failed to converge on %i particles.", count);
766
#endif
767

768
769
770
    /* Be clean */
    free(pid);
  }
771

772
  if (timer) TIMER_TOC(timer_do_ghost);
773
774
}

775
/**
776
 * @brief Unskip any tasks associated with active cells.
777
778
 *
 * @param c The cell.
779
 * @param e The engine.
780
 */
781
static void runner_do_unskip(struct cell *c, struct engine *e) {
782

783
784
785
  /* Ignore empty cells. */
  if (c->count == 0 && c->gcount == 0) return;

786
  /* Unskip any active tasks. */
787
  if (cell_is_active(c, e)) {
788
    const int forcerebuild = cell_unskip_tasks(c, &e->sched);
789
790
791
    if (forcerebuild) atomic_inc(&e->forcerebuild);
  }

792
  /* Recurse */
793
794
  if (c->split) {
    for (int k = 0; k < 8; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
795
      if (c->progeny[k] != NULL) {
Matthieu Schaller's avatar
Matthieu Schaller committed
796
        struct cell *cp = c->progeny[k];
797
        runner_do_unskip(cp, e);
798
799
800
      }
    }
  }
801
}
802

803
/**
804
 * @brief Mapper function to unskip active tasks.
805
806
807
808
809
 *
 * @param map_data An array of #cell%s.
 * @param num_elements Chunk size.
 * @param extra_data Pointer to an #engine.
 */
810
811
void runner_do_unskip_mapper(void *map_data, int num_elements,
                             void *extra_data) {
812

813
814
  struct engine *e = (struct engine *)extra_data;
  struct cell *cells = (struct cell *)map_data;
815

816
817
  for (int ind = 0; ind < num_elements; ind++) {
    struct cell *c = &cells[ind];
818
    if (c != NULL) runner_do_unskip(c, e);
819
  }
820
}
821
822
823
824
825
826
827
/**
 * @brief Drift particles in real space.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
828
void runner_do_drift_particles(struct runner *r, struct cell *c, int timer) {
829

830
  TIMER_TIC;
Matthieu Schaller's avatar
Matthieu Schaller committed
831

832
  cell_drift_particles(c, r->e);
833

834
  if (timer) TIMER_TOC(timer_drift);
835
}
836

837
838
839
840
841
842
843
/**
 * @brief Perform the first half-kick on all the active particles in a cell.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
844
void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
845

846
847
848
849
  const struct engine *e = r->e;
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
  struct gpart *restrict gparts = c->gparts;
850
  struct spart *restrict sparts = c->sparts;
851
852
  const int count = c->count;
  const int gcount = c->gcount;
853
  const int scount = c->scount;
854
  const integertime_t ti_current = e->ti_current;
855
  const double timeBase = e->timeBase;
856

857
858
859
  TIMER_TIC;

  /* Anything to do here? */
860
  if (!cell_is_starting(c, e)) return;
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_kick1(r, c->progeny[k], 0);
  } else {

    /* Loop over the parts in this cell. */
    for (int k = 0; k < count; k++) {

      /* Get a handle on the part. */
      struct part *restrict p = &parts[k];
      struct xpart *restrict xp = &xparts[k];

      /* If particle needs to be kicked */
876
      if (part_is_starting(p, e)) {
877
878
879

        const integertime_t ti_step = get_integer_timestep(p->time_bin);
        const integertime_t ti_begin =
880
            get_integer_time_begin(ti_current + 1, p->time_bin);
881

882
883
#ifdef SWIFT_DEBUG_CHECKS
        const integertime_t ti_end =
Matthieu Schaller's avatar
Matthieu Schaller committed
884
            get_integer_time_end(ti_current + 1, p->time_bin);
885

Matthieu Schaller's avatar
Matthieu Schaller committed
886
        if (ti_begin != ti_current)
887
888
889
890
          error(
              "Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
              "ti_step=%lld time_bin=%d ti_current=%lld",
              ti_end, ti_begin, ti_step, p->time_bin, ti_current);
891
892
#endif

893
        /* do the kick */
894
        kick_part(p, xp, ti_begin, ti_begin + ti_step / 2, timeBase);
895
896
897
      }
    }

898
    /* Loop over the gparts in this cell. */
899
900
901
902
903
904
    for (int k = 0; k < gcount; k++) {

      /* Get a handle on the part. */
      struct gpart *restrict gp = &gparts[k];

      /* If the g-particle has no counterpart and needs to be kicked */
905
      if (gp->type == swift_type_dark_matter && gpart_is_starting(gp, e)) {
906

907
908
        const integertime_t ti_step = get_integer_timestep(gp->time_bin);
        const integertime_t ti_begin =
909
            get_integer_time_begin(ti_current + 1, gp->time_bin);
910

911
912
#ifdef SWIFT_DEBUG_CHECKS
        const integertime_t ti_end =
Matthieu Schaller's avatar
Matthieu Schaller committed
913
            get_integer_time_end(ti_current + 1, gp->time_bin);
914

Matthieu Schaller's avatar
Matthieu Schaller committed
915
916
917
918
919
        if (ti_begin != ti_current)
          error(
              "Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
              "ti_step=%lld time_bin=%d ti_current=%lld",
              ti_end, ti_begin, ti_step, gp->time_bin, ti_current);
920
921
#endif

922
        /* do the kick */
923
        kick_gpart(gp, ti_begin, ti_begin + ti_step / 2, timeBase);
924
925
      }
    }
926
927
928
929
930
931
932
933

    /* Loop over the star particles in this cell. */
    for (int k = 0; k < scount; k++) {

      /* Get a handle on the s-part. */
      struct spart *restrict sp = &sparts[k];

      /* If particle needs to be kicked */
934
      if (spart_is_starting(sp, e)) {
935
936
937

        const integertime_t ti_step = get_integer_timestep(sp->time_bin);
        const integertime_t ti_begin =
938
            get_integer_time_begin(ti_current + 1, sp->time_bin);
939
940
941

#ifdef SWIFT_DEBUG_CHECKS
        const integertime_t ti_end =
942
            get_integer_time_end(ti_current + 1, sp->time_bin);
943