runner.c 38.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. */
Matthieu Schaller's avatar
Matthieu Schaller committed
41
#include "approx_math.h"
42
#include "atomic.h"
43
#include "cell.h"
44
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
45
#include "cooling.h"
46
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
47
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
48
#include "engine.h"
49
#include "error.h"
50
51
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
52
#include "hydro_properties.h"
53
#include "kick.h"
54
#include "minmax.h"
55
#include "potential.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
89
    {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},
    {0.0, 0.0, 1.0},
Matthieu Schaller's avatar
Matthieu Schaller committed
90
};
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
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
164
  const struct external_potential *potential = r->e->external_potential;
165
  const struct phys_const *constants = r->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
171
172
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

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 (gp->ti_end <= ti_current) {
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 int ti_current = r->e->ti_current;
497
498

  TIMER_TIC;
499

500
501
502
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

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 (p->ti_end <= ti_current) {
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
530
531

      if (gp->ti_end <= ti_current) {

        /* 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
553
554
  struct part *restrict parts = c->parts;
  const int count = c->count;
  const int ti_current = r->e->ti_current;

555
556
557
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
  /* 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];

      if (p->ti_end <= ti_current) {

        /* 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
  const int ti_current = r->e->ti_current;
  const double timeBase = r->e->timeBase;
598
599
600
601
602
603
604
  const float target_wcount = r->e->hydro_properties->target_neighbours;
  const float max_wcount =
      target_wcount + r->e->hydro_properties->delta_neighbours;
  const float min_wcount =
      target_wcount - r->e->hydro_properties->delta_neighbours;
  const int max_smoothing_iter =
      r->e->hydro_properties->max_smoothing_iterations;
605

606
607
  TIMER_TIC;

608
609
610
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

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 (p->ti_end <= ti_current) {
640

641
        /* Finish the density calculation */
642
        hydro_end_density(p, ti_current);
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
/**
751
 * @brief Drift particles and g-particles in a cell forward in time
752
753
 *
 * @param c The cell.
754
 * @param e The engine.
755
 */
756
static void runner_do_drift(struct cell *c, struct engine *e) {
757

758
  const double timeBase = e->timeBase;
759
  const int ti_old = c->ti_old;
760
  const int ti_current = e->ti_current;
761
762
763
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;
764
765
  const struct phys_const *constants = e->physical_constants;
  const struct external_potential *potential = e->external_potential;
766

767
768
769
  /* Do we need to drift ? */
  if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return;

770
771
772
  /* Check that we are actually going to move forward. */
  if (ti_current == ti_old) return;

773
774
775
776
  /* 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;
777
  double e_kin = 0.0, e_int = 0.0, e_pot_self = 0.0, e_pot_ext = 0.0, e_rad = 0.0;
778
  double entropy = 0.0, mass = 0.0;
779
780
781
  double mom[3] = {0.0, 0.0, 0.0};
  double ang_mom[3] = {0.0, 0.0, 0.0};

782
783
  /* No children? */
  if (!c->split) {
784

785
    /* Loop over all the g-particles in the cell */
786
    const size_t nr_gparts = c->gcount;
Peter W. Draper's avatar
Peter W. Draper committed
787
788
    for (size_t k = 0; k < nr_gparts; k++) {

789
      /* Get a handle on the gpart. */
790
      struct gpart *const gp = &gparts[k];
791
792

      /* Drift... */
793
      drift_gpart(gp, dt, timeBase, ti_old, ti_current);
794
795
796
797
798

      /* Compute (square of) motion since last cell construction */
      const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
                        gp->x_diff[1] * gp->x_diff[1] +
                        gp->x_diff[2] * gp->x_diff[2];
799
      dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
800
801
802
    }

    /* Loop over all the particles in the cell (more work for these !) */
Peter W. Draper's avatar
Peter W. Draper committed
803
804
    const size_t nr_parts = c->count;
    for (size_t k = 0; k < nr_parts; k++) {
805

806
      /* Get a handle on the part. */
807
808
      struct part *const p = &parts[k];
      struct xpart *const xp = &xparts[k];
809

810
      /* Drift... */
811
      drift_part(p, xp, dt, timeBase, ti_old, ti_current);
812

813
      /* Compute (square of) motion since last cell construction */
814
815
816
      const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
                        xp->x_diff[1] * xp->x_diff[1] +
                        xp->x_diff[2] * xp->x_diff[2];
817
      dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
Matthieu Schaller's avatar
Matthieu Schaller committed
818
819

      /* Maximal smoothing length */
820
      h_max = (h_max > p->h) ? h_max : p->h;
821
822
823
824
825
826
827
828
829

      /* Now collect quantities for statistics */

      const float half_dt =
          (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
      const double x[3] = {p->x[0], p->x[1], p->x[2]};
      const float v[3] = {xp->v_full[0] + p->a_hydro[0] * half_dt,
                          xp->v_full[1] + p->a_hydro[1] * half_dt,
                          xp->v_full[2] + p->a_hydro[2] * half_dt};
830
831

      const float m = hydro_get_mass(p);
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847

      /* Collect mass */
      mass += m;

      /* Collect momentum */
      mom[0] += m * v[0];
      mom[1] += m * v[1];
      mom[2] += m * v[2];

      /* Collect angular momentum */
      ang_mom[0] += m * (x[1] * v[2] - x[2] * v[1]);
      ang_mom[1] += m * (x[2] * v[0] - x[0] * v[2]);
      ang_mom[2] += m * (x[0] * v[1] - x[1] * v[0]);

      /* Collect energies. */
      e_kin += 0.5 * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
848
849
      e_pot_self += 0.;
      e_pot_ext += m * external_gravity_get_potential_energy(potential,constants,p);
850
      e_int += m * hydro_get_internal_energy(p, half_dt);
851
      e_rad += cooling_get_radiated_energy(xp);
852
853

      /* Collect entropy */
854
      entropy += m * hydro_get_entropy(p, half_dt);
855
    }
856
857

    /* Now, get the maximal particle motion from its square */
858
    dx_max = sqrtf(dx2_max);
859
  }
860

Matthieu Schaller's avatar
Matthieu Schaller committed
861
862
863
  /* Otherwise, aggregate data from children. */
  else {

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

869
870
        /* Recurse. */
        runner_do_drift(cp, e);
Matthieu Schaller's avatar
Matthieu Schaller committed
871

872
873
        dx_max = max(dx_max, cp->dx_max);
        h_max = max(h_max, cp->h_max);
874
875
876
        mass += cp->mass;
        e_kin += cp->e_kin;
        e_int += cp->e_int;
877
878
        e_pot_self += cp->e_pot_self;
	e_pot_ext += cp->e_pot_ext;
879
        e_rad += cp->e_rad;
880
        entropy += cp->entropy;
881
882
883
884
885
886
        mom[0] += cp->mom[0];
        mom[1] += cp->mom[1];
        mom[2] += cp->mom[2];
        ang_mom[0] += cp->ang_mom[0];
        ang_mom[1] += cp->ang_mom[1];
        ang_mom[2] += cp->ang_mom[2];
Matthieu Schaller's avatar
Matthieu Schaller committed
887
888
889
890
891
892
      }
  }

  /* Store the values */
  c->h_max = h_max;
  c->dx_max = dx_max;
893
894
895
  c->mass = mass;
  c->e_kin = e_kin;
  c->e_int = e_int;
896
897
  c->e_pot_self = e_pot_self;
  c->e_pot_ext = e_pot_ext;
898
  c->e_rad = e_rad;
899
  c->entropy = entropy;
900
901
902
903
904
905
  c->mom[0] = mom[0];
  c->mom[1] = mom[1];
  c->mom[2] = mom[2];
  c->ang_mom[0] = ang_mom[0];
  c->ang_mom[1] = ang_mom[1];
  c->ang_mom[2] = ang_mom[2];
906

907
908
  /* Update the time of the last drift */
  c->ti_old = ti_current;
909
}
910

911
912
913
914
915
916
917
918
919
920
/**
 * @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) {
921

922
923
924
925
926
927
928
  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];

    /* Only drift local particles. */
Peter W. Draper's avatar
Peter W. Draper committed
929
    if (c != NULL && c->nodeID == e->nodeID) runner_do_drift(c, e);
930
  }
931
}
932

933
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
934
935
 * @brief Kick particles in momentum space and collect statistics (fixed
 * time-step case)
936
937
938
 *
 * @param r The runner thread.
 * @param c The cell.
939
 * @param timer Are we timing this ?
940
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
941
void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
942

Matthieu Schaller's avatar
Matthieu Schaller committed
943
  const double global_dt = r->e->dt_max;
944
  const double timeBase = r->e->timeBase;
Matthieu Schaller's avatar
Matthieu Schaller committed
945
  const int count = c->count;
946
  const int gcount = c->gcount;
Matthieu Schaller's avatar
Matthieu Schaller committed
947
948
949
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
  struct gpart *restrict gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
950
  const double const_G = r->e->physical_constants->const_newton_G;
Matthieu Schaller's avatar
Matthieu Schaller committed
951

952
  int updated = 0, g_updated = 0;
953
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
954
955
956

  TIMER_TIC

Tom Theuns's avatar
Tom Theuns committed
957
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
958
  OUT;
Tom Theuns's avatar
Tom Theuns committed
959
960
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
961
962
963
  /* The new time-step */
  const int new_dti = global_dt / timeBase;

964
965
966
  /* No children? */
  if (!c->split) {

967
    /* Loop over the g-particles and kick everyone. */
968
969
970
    for (int k = 0; k < gcount; k++) {

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

973
      /* If the g-particle has no counterpart */
974
      if (gp->id_or_neg_offset > 0) {
975

Matthieu Schaller's avatar
Matthieu Schaller committed
976
        /* First, finish the force calculation */
977
        gravity_end_force(gp, const_G);
978

979
980
        /* Kick the g-particle forward */
        kick_gpart(gp, new_dti, timeBase);
981

Matthieu Schaller's avatar
Matthieu Schaller committed
982
983
        /* Number of updated g-particles */
        g_updated++;
984

Matthieu Schaller's avatar
Matthieu Schaller committed
985
986
987
988
989
        /* 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);
      }
    }
990

Matthieu Schaller's avatar
Matthieu Schaller committed
991
    /* Now do the hydro ones... */
992

993
    /* Loop over the particles and kick everyone. */
Matthieu Schaller's avatar
Matthieu Schaller committed
994
    for (int k = 0; k < count; k++) {
995

Matthieu Schaller's avatar
Matthieu Schaller committed
996
      /* Get a handle on the part. */
997
998
      struct part *restrict p = &parts[k];
      struct xpart *restrict xp = &xparts[k];
999

Matthieu Schaller's avatar
Matthieu Schaller committed