runner.c 35.7 KB
Newer Older
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
5
6
7
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
 *               2016 John A. Regan (john.a.regan@durham.ac.uk)
 *                    Tom Theuns (tom.theuns@durham.ac.uk)
8
 *
9
10
11
12
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
13
 *
14
15
16
17
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
18
 *
19
20
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
21
 *
22
 ******************************************************************************/
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"
56
#include "scheduler.h"
57
#include "sourceterms.h"
58
59
60
#include "space.h"
#include "task.h"
#include "timers.h"
61
#include "timestep.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
62

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

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

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

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

/* Does the axis need flipping ? */
93
94
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
95

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

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

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

112
/* Import the gravity loop functions. */
113
#include "runner_doiact_fft.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
114
#include "runner_doiact_grav.h"
115

Tom Theuns's avatar
Tom Theuns committed
116
/**
Tom Theuns's avatar
Tom Theuns committed
117
 * @brief Perform source terms
Tom Theuns's avatar
Tom Theuns committed
118
119
120
121
122
123
124
 *
 * @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;
125
  const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
Tom Theuns's avatar
Tom Theuns committed
126
  const double cell_width[3] = {c->width[0], c->width[1], c->width[2]};
Tom Theuns's avatar
Tom Theuns committed
127
  struct sourceterms *sourceterms = r->e->sourceterms;
128
  const int dimen = 3;
Tom Theuns's avatar
Tom Theuns committed
129
130
131
132
133
134
135
136
137
138

  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);
    return;
  }

Tom Theuns's avatar
Tom Theuns committed
139
140
141
142
143
144
145
  if (count > 0) {

    /* 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
146
147
    }
  }
Tom Theuns's avatar
Tom Theuns committed
148
149
150
151

  if (timer) TIMER_TOC(timer_dosource);
}

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

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

168
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
169

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

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

#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
181
  OUT;
Tom Theuns's avatar
Tom Theuns committed
182
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
183

184
  /* Loop over the gparts in this cell. */
Matthieu Schaller's avatar
Matthieu Schaller committed
185
  for (int i = 0; i < gcount; i++) {
186

187
    /* Get a direct pointer on the part. */
188
    struct gpart *restrict gp = &gparts[i];
189
190

    /* Is this part within the time step? */
191
    if (gpart_is_active(gp, e)) {
Matthieu Schaller's avatar
Matthieu Schaller committed
192

193
      external_gravity_acceleration(time, potential, constants, gp);
194
    }
195
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
196

197
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
198
199
}

Stefan Arridge's avatar
Stefan Arridge committed
200
/**
201
202
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
203
204
205
206
207
208
209
210
 *
 * @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;
211
  struct xpart *restrict xparts = c->xparts;
Stefan Arridge's avatar
Stefan Arridge committed
212
213
  const int count = c->count;
  const int ti_current = r->e->ti_current;
214
  const struct cooling_function_data *cooling_func = r->e->cooling_func;
Stefan Arridge's avatar
Stefan Arridge committed
215
  const struct phys_const *constants = r->e->physical_constants;
216
  const struct UnitSystem *us = r->e->internalUnits;
Stefan Arridge's avatar
Stefan Arridge committed
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
  const double timeBase = r->e->timeBase;

  TIMER_TIC;

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

#ifdef TASK_VERBOSE
  OUT;
#endif

Stefan Arridge's avatar
Stefan Arridge committed
232
  /* Loop over the parts in this cell. */
Stefan Arridge's avatar
Stefan Arridge committed
233
234
235
236
  for (int i = 0; i < count; i++) {

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

Stefan Arridge's avatar
Stefan Arridge committed
239
240
    /* Kick has already updated ti_end, so need to check ti_begin */
    if (p->ti_begin == ti_current) {
241
242
243

      const double dt = (p->ti_end - p->ti_begin) * timeBase;

244
      cooling_cool_part(constants, us, cooling_func, p, xp, dt);
Stefan Arridge's avatar
Stefan Arridge committed
245
246
247
248
249
250
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
251
252
253
254
255
256
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
257
void runner_do_sort_ascending(struct entry *sort, int N) {
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
307
308
309
310
311

  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
312
        }
313
314
315
316
317
318
319
320
321
322
323
324
      } 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
325
    }
326
327
328
  }
}

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

  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
346
347
  float buff[8];
  double px[3];
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371

  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;
372
      if (missing) runner_do_sort(r, c->progeny[k], missing, 0);
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
420
421
422
423
424
    }

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

427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
      } /* 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
451
452
453
          sort[j * (count + 1) + k].d = px[0] * runner_shift[j][0] +
                                        px[1] * runner_shift[j][1] +
                                        px[2] * runner_shift[j][2];
454
        }
455
    }
456
457
458
459
460
461

    /* 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;
462
        runner_do_sort_ascending(&sort[j * (count + 1)], count);
463
464
465
466
        c->sorted |= (1 << j);
      }
  }

467
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
468
  /* Verify the sorting. */
469
470
471
472
473
474
475
476
477
478
  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
479
480
481
482

  if (clock) TIMER_TOC(timer_dosort);
}

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

Matthieu Schaller's avatar
Matthieu Schaller committed
492
493
  struct part *restrict parts = c->parts;
  struct gpart *restrict gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
494
  const int count = c->count;
495
  const int gcount = c->gcount;
496
  const struct engine *e = r->e;
497
498

  TIMER_TIC;
499

500
  /* Anything to do here? */
501
  if (!cell_is_active(c, e)) return;
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
    return;
508
509
  } else {

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

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

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

518
519
        /* Get ready for a density calculation */
        hydro_init_part(p);
520
      }
521
    }
522
523
524
525
526

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

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

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

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

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

540
/**
541
542
543
544
545
546
547
 * @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.
 */
void runner_do_extra_ghost(struct runner *r, struct cell *c) {
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
  /* Anything to do here? */
556
  if (!cell_is_active(c, e)) return;
557

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

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

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

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

584
/**
585
586
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
587
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
588
 * @param r The runner thread.
589
 * @param c The cell.
590
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
591
void runner_do_ghost(struct runner *r, struct cell *c) {
592

Matthieu Schaller's avatar
Matthieu Schaller committed
593
594
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
595
  int redo, count = c->count;
596
597
598
599
  const struct engine *e = r->e;
  const int ti_current = e->ti_current;
  const double timeBase = e->timeBase;
  const float target_wcount = e->hydro_properties->target_neighbours;
600
  const float max_wcount =
601
      target_wcount + e->hydro_properties->delta_neighbours;
602
  const float min_wcount =
603
604
      target_wcount - e->hydro_properties->delta_neighbours;
  const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
605

606
607
  TIMER_TIC;

608
  /* Anything to do here? */
609
  if (!cell_is_active(c, e)) return;
610

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

  /* Init the IDs that have to be updated. */
619
620
621
  int *pid = NULL;
  if ((pid = malloc(sizeof(int) * count)) == NULL)
    error("Can't allocate memory for pid.");
Matthieu Schaller's avatar
Matthieu Schaller committed
622
  for (int k = 0; k < count; k++) pid[k] = k;
623
624

  /* While there are particles that need to be updated... */
625
  for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
Matthieu Schaller's avatar
Matthieu Schaller committed
626
       num_reruns++) {
627
628
629

    /* Reset the redo-count. */
    redo = 0;
630

631
    /* Loop over the parts in this cell. */
Matthieu Schaller's avatar
Matthieu Schaller committed
632
    for (int i = 0; i < count; i++) {
633
634

      /* Get a direct pointer on the part. */
635
636
      struct part *restrict p = &parts[pid[i]];
      struct xpart *restrict xp = &xparts[pid[i]];
637
638

      /* Is this part within the timestep? */
639
      if (part_is_active(p, e)) {
640

641
        /* Finish the density calculation */
642
        hydro_end_density(p);
643

Matthieu Schaller's avatar
Matthieu Schaller committed
644
645
        float h_corr = 0.f;

646
        /* If no derivative, double the smoothing length. */
647
        if (p->density.wcount_dh == 0.0f) h_corr = p->h;
648
649
650

        /* Otherwise, compute the smoothing length update (Newton step). */
        else {
651
          h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
652
653

          /* Truncate to the range [ -p->h/2 , p->h ]. */
654
655
          h_corr = (h_corr < p->h) ? h_corr : p->h;
          h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
Pedro Gonnet's avatar
Pedro Gonnet committed
656
        }
657
658

        /* Did we get the right number density? */
659
        if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
660
661
662
663

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

664
          /* Flag for another round of fun */
665
666
          pid[redo] = pid[i];
          redo += 1;
667

668
669
          /* Re-initialise everything */
          hydro_init_part(p);
670

671
          /* Off we go ! */
672
          continue;
673
674
        }

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

677
        /* As of here, particle force variables will be set. */
678

679
        /* Compute variables required for the force loop */
680
        hydro_prepare_force(p, xp, ti_current, timeBase);
681

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

685
686
        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
687
688
689
      }
    }

690
691
692
    /* We now need to treat the particles whose smoothing length had not
     * converged again */

693
694
695
696
697
    /* Re-set the counter for the next loop (potentially). */
    count = redo;
    if (count > 0) {

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

700
701
        /* 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
702

703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
          /* Self-interaction? */
          if (l->t->type == task_type_self)
            runner_doself_subset_density(r, finger, parts, pid, count);

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

            /* 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);

          }

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

          /* Otherwise, sub-pair interaction? */
          else if (l->t->type == task_type_sub_pair) {
727
728
729
730
731
732
733
734
735
736
737

            /* 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);
          }
        }
      }
738
    }
739
  }
740

Matthieu Schaller's avatar
Matthieu Schaller committed
741
742
  if (count)
    message("Smoothing length failed to converge on %i particles.", count);
743

744
745
746
  /* Be clean */
  free(pid);

Matthieu Schaller's avatar
Matthieu Schaller committed
747
  TIMER_TOC(timer_do_ghost);
748
749
}

750
/**
Peter W. Draper's avatar
Peter W. Draper committed
751
752
 * @brief Drift particles and g-particles in a cell forward in time,
 *              unskipping any tasks associated with active cells.
753
754
 *
 * @param c The cell.
755
 * @param e The engine.
Peter W. Draper's avatar
Peter W. Draper committed
756
757
 * @param drift whether to actually drift the particles, will not be
 *              necessary for non-local cells.
758
 */
759
static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
760

761
  /* Unskip any active tasks. */
762
  if (cell_is_active(c, e)) {
763
    const int forcerebuild = cell_unskip_tasks(c, &e->sched);
764
765
766
767
768
    if (forcerebuild) atomic_inc(&e->forcerebuild);
  }

  /* Do we really need to drift? */
  if (drift) {
769
    if (!e->drift_all && !cell_is_drift_needed(c, e)) return;
770
771
  } else {

772
    /* Not drifting, but may still need to recurse for task un-skipping. */
773
774
775
776
777
778
779
780
781
782
783
    if (c->split) {
      for (int k = 0; k < 8; k++) {
        if (c->progeny[k] != NULL) {
          struct cell *cp = c->progeny[k];
          runner_do_drift(cp, e, 0);
        }
      }
    }
    return;
  }

784
785
786
  /* Now, we can drift */

  /* Get some information first */
787
788
  const double timeBase = e->timeBase;
  const int ti_old = c->ti_old;
789
  const int ti_current = e->ti_current;
790
791
792
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;
793
794
795
796

  /* Drift from the last time the cell was drifted to the current time */
  const double dt = (ti_current - ti_old) * timeBase;
  float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
797

798
799
  /* No children? */
  if (!c->split) {
800

801
    /* Check that we are actually going to move forward. */
802
    if (ti_current > ti_old) {
Peter W. Draper's avatar
Peter W. Draper committed
803

804
805
806
      /* Loop over all the g-particles in the cell */
      const size_t nr_gparts = c->gcount;
      for (size_t k = 0; k < nr_gparts; k++) {
807

808
809
        /* Get a handle on the gpart. */
        struct gpart *const gp = &gparts[k];
810

811
812
        /* Drift... */
        drift_gpart(gp, dt, timeBase, ti_old, ti_current);
813

814
815
        /* Compute (square of) motion since last cell construction */
        const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
816
817
                          gp->x_diff[1] * gp->x_diff[1] +
                          gp->x_diff[2] * gp->x_diff[2];
818
819
        dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
      }
820

821
      /* Loop over all the particles in the cell */
822
823
      const size_t nr_parts = c->count;
      for (size_t k = 0; k < nr_parts; k++) {
824

825
826
827
        /* Get a handle on the part. */
        struct part *const p = &parts[k];
        struct xpart *const xp = &xparts[k];
828

829
830
        /* Drift... */
        drift_part(p, xp, dt, timeBase, ti_old, ti_current);
Matthieu Schaller's avatar
Matthieu Schaller committed
831

832
833
        /* Compute (square of) motion since last cell construction */
        const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
834
835
                          xp->x_diff[1] * xp->x_diff[1] +
                          xp->x_diff[2] * xp->x_diff[2];
836
        dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
837

838
839
840
        /* Maximal smoothing length */
        h_max = (h_max > p->h) ? h_max : p->h;
      }
841

842
843
      /* Now, get the maximal particle motion from its square */
      dx_max = sqrtf(dx2_max);
844

845
    } /* Check that we are actually going to move forward. */
846
847
848
849
850
851

    else {
      /* ti_old == ti_current, just keep the current cell values. */
      h_max = c->h_max;
      dx_max = c->dx_max;
    }
852
  }
853

Matthieu Schaller's avatar
Matthieu Schaller committed
854
855
856
  /* Otherwise, aggregate data from children. */
  else {

857
    /* Loop over the progeny and collect their data. */
Matthieu Schaller's avatar
Matthieu Schaller committed
858
859
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
860
        struct cell *cp = c->progeny[k];
861

862
        /* Recurse. */
863
        runner_do_drift(cp, e, drift);
864
865
        dx_max = max(dx_max, cp->dx_max);
        h_max = max(h_max, cp->h_max);
Matthieu Schaller's avatar
Matthieu Schaller committed
866
867
868
869
      }
  }

  /* Store the values */
870
871
872
873
874
  c->h_max = h_max;
  c->dx_max = dx_max;

  /* Update the time of the last drift */
  c->ti_old = ti_current;
875
}
876

877
878
879
880
881
882
883
884
885
886
/**
 * @brief Mapper function to drift particles and g-particles forward in time.
 *
 * @param map_data An array of #cell%s.
 * @param num_elements Chunk size.
 * @param extra_data Pointer to an #engine.
 */

void runner_do_drift_mapper(void *map_data, int num_elements,
                            void *extra_data) {
887

888
889
890
891
892
  struct engine *e = (struct engine *)extra_data;
  struct cell *cells = (struct cell *)map_data;

  for (int ind = 0; ind < num_elements; ind++) {
    struct cell *c = &cells[ind];
893
#ifdef WITH_MPI
894
    if (c != NULL) runner_do_drift(c, e, (c->nodeID == e->nodeID));
895
896
897
#else
    if (c != NULL) runner_do_drift(c, e, 1);
#endif
898
  }
899
}
900

Matthieu Schaller's avatar
Matthieu Schaller committed
901
902
903
/**
 * @brief Kick particles in momentum space and collect statistics (floating
 * time-step case)
904
905
906
 *
 * @param r The runner thread.
 * @param c The cell.
907
 * @param timer Are we timing this ?
908
 */
909
void runner_do_kick(struct runner *r, struct cell *c, int timer) {
910

911
912
  const struct engine *e = r->e;
  const double timeBase = e->timeBase;
Matthieu Schaller's avatar
Matthieu Schaller committed
913
  const int count = c->count;
914
  const int gcount = c->gcount;
Matthieu Schaller's avatar
Matthieu Schaller committed
915
916
917
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
  struct gpart *restrict gparts = c->gparts;
918
  const double const_G = e->physical_constants->const_newton_G;
Matthieu Schaller's avatar
Matthieu Schaller committed
919

920
  TIMER_TIC;
921

922
  /* Anything to do here? */
923
  if (!cell_is_active(c, e)) {
924
925
926
927
    c->updated = 0;
    c->g_updated = 0;
    return;
  }
928

Tom Theuns's avatar
Tom Theuns committed
929
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
930
  OUT;
Tom Theuns's avatar
Tom Theuns committed
931
932
#endif

933
934
935
  int updated = 0, g_updated = 0;
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;

936
937
938
  /* No children? */
  if (!c->split) {

939
940
941
942
    /* Loop over the g-particles and kick the active ones. */
    for (int k = 0; k < gcount; k++) {

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

      /* If the g-particle has no counterpart and needs to be kicked */
946
      if (gp->id_or_neg_offset > 0) {
947

948
        if (gpart_is_active(gp, e)) {
949

Matthieu Schaller's avatar
Matthieu Schaller committed
950
          /* First, finish the force calculation */
Matthieu Schaller's avatar
Matthieu Schaller committed
951
          gravity_end_force(gp, const_G);
952

953
          /* Compute the next timestep */
954
          const int new_dti = get_gpart_timestep(gp, e);
955

Matthieu Schaller's avatar
Matthieu Schaller committed
956
          /* Now we have a time step, proceed with the kick */
957
          kick_gpart(gp, new_dti, timeBase);
Matthieu Schaller's avatar
Matthieu Schaller committed
958
959
960
961
962
963
964
965

          /* Number of updated g-particles */
          g_updated++;
        }

        /* Minimal time for next end of time-step */
        ti_end_min = min(gp->ti_end, ti_end_min);
        ti_end_max = max(gp->ti_end, ti_end_max);
966
967
968
969
970
      }
    }

    /* Now do the hydro ones... */

971
    /* Loop over the particles and kick the active ones. */
Matthieu Schaller's avatar
Matthieu Schaller committed
972
    for (int k = 0; k < count; k++) {
973
974

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

      /* If particle needs to be kicked */
979
      if (part_is_active(p, e)) {
980
981

        /* First, finish the force loop */
982
        hydro_end_force(p);
Matthieu Schaller's avatar
Matthieu Schaller committed
983
        if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
984