runner.c 32.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. */
Matthieu Schaller's avatar
Matthieu Schaller committed
41
#include "approx_math.h"
42
#include "atomic.h"
43
#include "cell.h"
44
#include "const.h"
45
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
46
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
47
#include "engine.h"
48
#include "error.h"
49
50
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
51
#include "hydro_properties.h"
52
#include "kick.h"
53
#include "minmax.h"
54
55
56
57
#include "scheduler.h"
#include "space.h"
#include "task.h"
#include "timers.h"
58
#include "timestep.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
59

60
/* Orientation of the cell pairs */
Matthieu Schaller's avatar
Matthieu Schaller committed
61
62
const double runner_shift[13][3] = {
    {5.773502691896258e-01, 5.773502691896258e-01, 5.773502691896258e-01},
63
    {7.071067811865475e-01, 7.071067811865475e-01, 0.0},
Matthieu Schaller's avatar
Matthieu Schaller committed
64
65
66
67
68
69
70
71
72
73
74
    {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
75
};
76
77

/* Does the axis need flipping ? */
78
79
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
80

81
/* Import the density loop functions. */
82
83
84
#define FUNCTION density
#include "runner_doiact.h"

85
/* Import the force loop functions. */
86
87
88
89
#undef FUNCTION
#define FUNCTION force
#include "runner_doiact.h"

90
/* Import the gravity loop functions. */
91
#include "runner_doiact_fft.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
92
#include "runner_doiact_grav.h"
93

Tom Theuns's avatar
Tom Theuns committed
94
95
96
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
97
98
99
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
100
 */
101
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
102

Matthieu Schaller's avatar
Matthieu Schaller committed
103
104
  struct gpart *restrict gparts = c->gparts;
  const int gcount = c->gcount;
105
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
106
  const struct external_potential *potential = r->e->external_potential;
107
  const struct phys_const *constants = r->e->physical_constants;
Matthieu Schaller's avatar
Matthieu Schaller committed
108

109
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
110
111
112

  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
113
    for (int k = 0; k < 8; k++)
114
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
Tom Theuns's avatar
Tom Theuns committed
115
116
117
118
    return;
  }

#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
119
  OUT;
Tom Theuns's avatar
Tom Theuns committed
120
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
121

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

125
    /* Get a direct pointer on the part. */
Matthieu Schaller's avatar
Matthieu Schaller committed
126
    struct gpart *restrict g = &gparts[i];
127
128
129

    /* Is this part within the time step? */
    if (g->ti_end <= ti_current) {
Matthieu Schaller's avatar
Matthieu Schaller committed
130

131
      external_gravity(potential, constants, g);
132
    }
133
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
134

135
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
136
137
}

Pedro Gonnet's avatar
Pedro Gonnet committed
138
139
140
141
142
143
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
144
void runner_do_sort_ascending(struct entry *sort, int N) {
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198

  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
199
        }
200
201
202
203
204
205
206
207
208
209
210
211
      } 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
212
    }
213
214
215
  }
}

Pedro Gonnet's avatar
Pedro Gonnet committed
216
217
218
219
220
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
221
 * @param flags Cell flag.
222
223
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
224
 */
225
void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
226
227
228
229
230
231
232

  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
233
234
  float buff[8];
  double px[3];
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258

  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;
259
      if (missing) runner_do_sort(r, c->progeny[k], missing, 0);
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
    }

    /* 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
312
        }
313

314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
      } /* 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
338
339
340
          sort[j * (count + 1) + k].d = px[0] * runner_shift[j][0] +
                                        px[1] * runner_shift[j][1] +
                                        px[2] * runner_shift[j][2];
341
        }
342
    }
343
344
345
346
347
348

    /* 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;
349
        runner_do_sort_ascending(&sort[j * (count + 1)], count);
350
351
352
353
        c->sorted |= (1 << j);
      }
  }

354
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
355
  /* Verify the sorting. */
356
357
358
359
360
361
362
363
364
365
  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
366
367
368
369

  if (clock) TIMER_TOC(timer_dosort);
}

370
371
372
373
374
/**
 * @brief Initialize the particles before the density calculation
 *
 * @param r The runner thread.
 * @param c The cell.
Matthieu Schaller's avatar
Matthieu Schaller committed
375
 * @param timer 1 if the time is to be recorded.
376
 */
377
void runner_do_init(struct runner *r, struct cell *c, int timer) {
378

Matthieu Schaller's avatar
Matthieu Schaller committed
379
380
  struct part *restrict parts = c->parts;
  struct gpart *restrict gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
381
  const int count = c->count;
382
  const int gcount = c->gcount;
383
  const int ti_current = r->e->ti_current;
384
385

  TIMER_TIC;
386

387
388
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
389
    for (int k = 0; k < 8; k++)
390
      if (c->progeny[k] != NULL) runner_do_init(r, c->progeny[k], 0);
391
    return;
392
393
  } else {

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

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

400
      if (p->ti_end <= ti_current) {
Matthieu Schaller's avatar
Matthieu Schaller committed
401

402
403
        /* Get ready for a density calculation */
        hydro_init_part(p);
404
      }
405
    }
406
407
408
409
410

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

      /* Get a direct pointer on the part. */
411
      struct gpart *restrict gp = &gparts[i];
412
413
414
415

      if (gp->ti_end <= ti_current) {

        /* Get ready for a density calculation */
416
        gravity_init_gpart(gp);
417
418
      }
    }
419
  }
420

Peter W. Draper's avatar
Peter W. Draper committed
421
  if (timer) TIMER_TOC(timer_init);
422
423
}

424
425
426
/**
 * @brief Intermediate task between density and force
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
427
 * @param r The runner thread.
428
 * @param c The cell.
429
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
430
void runner_do_ghost(struct runner *r, struct cell *c) {
431

Matthieu Schaller's avatar
Matthieu Schaller committed
432
433
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
434
  int redo, count = c->count;
435
436
  const int ti_current = r->e->ti_current;
  const double timeBase = r->e->timeBase;
437
438
439
440
441
442
443
  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;
444

445
446
  TIMER_TIC;

447
448
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
449
    for (int k = 0; k < 8; k++)
Matthieu Schaller's avatar
Matthieu Schaller committed
450
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k]);
451
452
453
454
    return;
  }

  /* Init the IDs that have to be updated. */
455
456
457
  int *pid = NULL;
  if ((pid = malloc(sizeof(int) * count)) == NULL)
    error("Can't allocate memory for pid.");
Matthieu Schaller's avatar
Matthieu Schaller committed
458
  for (int k = 0; k < count; k++) pid[k] = k;
459
460

  /* While there are particles that need to be updated... */
461
  for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
Matthieu Schaller's avatar
Matthieu Schaller committed
462
       num_reruns++) {
463
464
465

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

467
    /* Loop over the parts in this cell. */
Matthieu Schaller's avatar
Matthieu Schaller committed
468
    for (int i = 0; i < count; i++) {
469
470

      /* Get a direct pointer on the part. */
471
472
      struct part *restrict p = &parts[pid[i]];
      struct xpart *restrict xp = &xparts[pid[i]];
473
474

      /* Is this part within the timestep? */
475
      if (p->ti_end <= ti_current) {
476

477
        /* Finish the density calculation */
478
        hydro_end_density(p, ti_current);
479

Matthieu Schaller's avatar
Matthieu Schaller committed
480
481
        float h_corr = 0.f;

482
        /* If no derivative, double the smoothing length. */
483
        if (p->density.wcount_dh == 0.0f) h_corr = p->h;
484
485
486

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

          /* Truncate to the range [ -p->h/2 , p->h ]. */
490
491
          h_corr = fminf(h_corr, p->h);
          h_corr = fmaxf(h_corr, -p->h * 0.5f);
Pedro Gonnet's avatar
Pedro Gonnet committed
492
        }
493
494

        /* Did we get the right number density? */
495
        if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
496
497
498
499

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

500
          /* Flag for another round of fun */
501
502
          pid[redo] = pid[i];
          redo += 1;
503

504
505
          /* Re-initialise everything */
          hydro_init_part(p);
506

507
          /* Off we go ! */
508
          continue;
509
510
        }

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

513
        /* As of here, particle force variables will be set. */
514

515
        /* Compute variables required for the force loop */
516
        hydro_prepare_force(p, xp, ti_current, timeBase);
517

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

521
522
        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
523
524
525
      }
    }

526
527
528
    /* We now need to treat the particles whose smoothing length had not
     * converged again */

529
530
531
532
533
    /* 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
534
      for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
Matthieu Schaller's avatar
Matthieu Schaller committed
535

536
537
        /* 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
538

539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
          /* 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);

          }

556
557
558
559
560
561
562
          /* 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) {
563
564
565
566
567
568
569
570
571
572
573

            /* 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);
          }
        }
      }
574
    }
575
  }
576

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

580
581
582
  /* Be clean */
  free(pid);

Matthieu Schaller's avatar
Matthieu Schaller committed
583
  TIMER_TOC(timer_do_ghost);
584
585
}

586
/**
587
 * @brief Drift particles and g-particles forward in time
588
589
590
 *
 * @param r The runner thread.
 * @param c The cell.
591
 * @param timer Are we timing this ?
592
 */
593
void runner_do_drift(struct runner *r, struct cell *c, int timer) {
594

595
596
  const double timeBase = r->e->timeBase;
  const double dt = (r->e->ti_current - r->e->ti_old) * timeBase;
597
598
  const int ti_old = r->e->ti_old;
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
599
600
601
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
  struct gpart *restrict gparts = c->gparts;
602
  float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
603

604
  double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
605
606
607
  double mom[3] = {0.0, 0.0, 0.0};
  double ang_mom[3] = {0.0, 0.0, 0.0};

608
  TIMER_TIC
609

Tom Theuns's avatar
Tom Theuns committed
610
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
611
  OUT;
Tom Theuns's avatar
Tom Theuns committed
612
#endif
613

614
615
  /* No children? */
  if (!c->split) {
616

617
    /* Loop over all the g-particles in the cell */
Peter W. Draper's avatar
Peter W. Draper committed
618
619
620
    const int nr_gparts = c->gcount;
    for (size_t k = 0; k < nr_gparts; k++) {

621
      /* Get a handle on the gpart. */
622
      struct gpart *restrict gp = &gparts[k];
623
624

      /* Drift... */
625
      drift_gpart(gp, dt, timeBase, ti_old, ti_current);
626
627
628
629
630
631

      /* 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];
      dx2_max = fmaxf(dx2_max, dx2);
632
633
634
    }

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

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

642
      /* Drift... */
643
      drift_part(p, xp, dt, timeBase, ti_old, ti_current);
644

645
      /* Compute (square of) motion since last cell construction */
646
647
648
      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];
649
      dx2_max = fmaxf(dx2_max, dx2);
Matthieu Schaller's avatar
Matthieu Schaller committed
650
651
652

      /* Maximal smoothing length */
      h_max = fmaxf(p->h, h_max);
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679

      /* 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};
      const float m = p->mass;

      /* 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]);
      e_pot += 0.;
680
      e_int += m * hydro_get_internal_energy(p, half_dt);
681
682

      /* Collect entropy */
683
      entropy += m * hydro_get_entropy(p, half_dt);
684
    }
685
686

    /* Now, get the maximal particle motion from its square */
687
    dx_max = sqrtf(dx2_max);
688
  }
689

Matthieu Schaller's avatar
Matthieu Schaller committed
690
691
692
693
694
695
  /* Otherwise, aggregate data from children. */
  else {

    /* Loop over the progeny. */
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
696
697

        /* Recurse */
698
        struct cell *restrict cp = c->progeny[k];
699
        runner_do_drift(r, cp, 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
700

701
        /* Collect */
702
703
        dx_max = fmaxf(dx_max, cp->dx_max);
        h_max = fmaxf(h_max, cp->h_max);
704
705
706
707
        mass += cp->mass;
        e_kin += cp->e_kin;
        e_int += cp->e_int;
        e_pot += cp->e_pot;
708
        entropy += cp->entropy;
709
710
711
712
713
714
        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
715
716
717
718
719
720
      }
  }

  /* Store the values */
  c->h_max = h_max;
  c->dx_max = dx_max;
721
722
723
724
  c->mass = mass;
  c->e_kin = e_kin;
  c->e_int = e_int;
  c->e_pot = e_pot;
725
  c->entropy = entropy;
726
727
728
729
730
731
  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];
732

Peter W. Draper's avatar
Peter W. Draper committed
733
  if (timer) TIMER_TOC(timer_drift);
734
}
735

736
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
737
738
 * @brief Kick particles in momentum space and collect statistics (fixed
 * time-step case)
739
740
741
 *
 * @param r The runner thread.
 * @param c The cell.
742
 * @param timer Are we timing this ?
743
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
744
void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
745

Matthieu Schaller's avatar
Matthieu Schaller committed
746
  const double global_dt = r->e->dt_max;
747
  const double timeBase = r->e->timeBase;
Matthieu Schaller's avatar
Matthieu Schaller committed
748
  const int count = c->count;
749
  const int gcount = c->gcount;
Matthieu Schaller's avatar
Matthieu Schaller committed
750
751
752
  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
753
  const double const_G = r->e->physical_constants->const_newton_G;
Matthieu Schaller's avatar
Matthieu Schaller committed
754

755
  int updated = 0, g_updated = 0;
756
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
757
758
759

  TIMER_TIC

Tom Theuns's avatar
Tom Theuns committed
760
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
761
  OUT;
Tom Theuns's avatar
Tom Theuns committed
762
763
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
764
765
766
  /* The new time-step */
  const int new_dti = global_dt / timeBase;

767
768
769
  /* No children? */
  if (!c->split) {

770
    /* Loop over the g-particles and kick everyone. */
771
772
773
    for (int k = 0; k < gcount; k++) {

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

776
      /* If the g-particle has no counterpart */
777
      if (gp->id_or_neg_offset > 0) {
778

Matthieu Schaller's avatar
Matthieu Schaller committed
779
        /* First, finish the force calculation */
780
        gravity_end_force(gp, const_G);
781

782
783
        /* Kick the g-particle forward */
        kick_gpart(gp, new_dti, timeBase);
784

Matthieu Schaller's avatar
Matthieu Schaller committed
785
786
        /* Number of updated g-particles */
        g_updated++;
787

Matthieu Schaller's avatar
Matthieu Schaller committed
788
789
790
791
792
        /* 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);
      }
    }
793

Matthieu Schaller's avatar
Matthieu Schaller committed
794
    /* Now do the hydro ones... */
795

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

Matthieu Schaller's avatar
Matthieu Schaller committed
799
      /* Get a handle on the part. */
800
801
      struct part *restrict p = &parts[k];
      struct xpart *restrict xp = &xparts[k];
802

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

807
808
      /* Kick the particle forward */
      kick_part(p, xp, new_dti, timeBase);
809

Matthieu Schaller's avatar
Matthieu Schaller committed
810
811
812
      /* Number of updated particles */
      updated++;
      if (p->gpart != NULL) g_updated++;
813

Matthieu Schaller's avatar
Matthieu Schaller committed
814
815
816
      /* Minimal time for next end of time-step */
      ti_end_min = min(p->ti_end, ti_end_min);
      ti_end_max = max(p->ti_end, ti_end_max);
Matthieu Schaller's avatar
Matthieu Schaller committed
817
818
    }
  }
819

Matthieu Schaller's avatar
Matthieu Schaller committed
820
821
  /* Otherwise, aggregate data from children. */
  else {
822

Matthieu Schaller's avatar
Matthieu Schaller committed
823
824
825
    /* Loop over the progeny. */
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
826
        struct cell *restrict cp = c->progeny[k];
827

Matthieu Schaller's avatar
Matthieu Schaller committed
828
829
        /* Recurse */
        runner_do_kick_fixdt(r, cp, 0);
830

Matthieu Schaller's avatar
Matthieu Schaller committed
831
832
833
834
835
836
837
        /* And aggregate */
        updated += cp->updated;
        g_updated += cp->g_updated;
        ti_end_min = min(cp->ti_end_min, ti_end_min);
        ti_end_max = max(cp->ti_end_max, ti_end_max);
      }
  }
838

Matthieu Schaller's avatar
Matthieu Schaller committed
839
840
841
842
843
  /* Store the values. */
  c->updated = updated;
  c->g_updated = g_updated;
  c->ti_end_min = ti_end_min;
  c->ti_end_max = ti_end_max;
844

Matthieu Schaller's avatar
Matthieu Schaller committed
845
846
  if (timer) TIMER_TOC(timer_kick);
}
847

Matthieu Schaller's avatar
Matthieu Schaller committed
848
849
850
/**
 * @brief Kick particles in momentum space and collect statistics (floating
 * time-step case)
851
852
853
 *
 * @param r The runner thread.
 * @param c The cell.
854
 * @param timer Are we timing this ?
855
 */
856
void runner_do_kick(struct runner *r, struct cell *c, int timer) {
857

858
859
  const struct engine *e = r->e;
  const double timeBase = e->timeBase;
860
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
861
  const int count = c->count;
862
  const int gcount = c->gcount;
Matthieu Schaller's avatar
Matthieu Schaller committed
863
864
865
  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
866
  const double const_G = r->e->physical_constants->const_newton_G;
Matthieu Schaller's avatar
Matthieu Schaller committed
867

868
  int updated = 0, g_updated = 0;
869
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
870
871

  TIMER_TIC
872

Tom Theuns's avatar
Tom Theuns committed
873
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
874
  OUT;
Tom Theuns's avatar
Tom Theuns committed
875
876
#endif

877
878
879
  /* No children? */
  if (!c->split) {

880
881
882
883
    /* Loop over the g-particles and kick the active ones. */
    for (int k = 0; k < gcount; k++) {

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
889
        if (gp->ti_end <= ti_current) {
890

Matthieu Schaller's avatar
Matthieu Schaller committed
891
          /* First, finish the force calculation */
Matthieu Schaller's avatar
Matthieu Schaller committed
892
          gravity_end_force(gp, const_G);
893

894
          /* Compute the next timestep */
895
          const int new_dti = get_gpart_timestep(gp, e);
896

Matthieu Schaller's avatar
Matthieu Schaller committed
897
          /* Now we have a time step, proceed with the kick */
898
          kick_gpart(gp, new_dti, timeBase);
Matthieu Schaller's avatar
Matthieu Schaller committed
899
900
901
902
903
904
905
906

          /* 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);
907
908
909
910
911
      }
    }

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

912
    /* Loop over the particles and kick the active ones. */
Matthieu Schaller's avatar
Matthieu Schaller committed
913
    for (int k = 0; k < count; k++) {
914
915

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

      /* If particle needs to be kicked */
Matthieu Schaller's avatar
Matthieu Schaller committed
920
      if (p->ti_end <= ti_current) {
921
922

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

Matthieu Schaller's avatar
Matthieu Schaller committed
926
        /* Compute the next timestep (hydro condition) */
927
        const int new_dti = get_part_timestep(p, xp, e);
928

Matthieu Schaller's avatar
Matthieu Schaller committed
929
        /* Now we have a time step, proceed with the kick */
930
        kick_part(p, xp, new_dti, timeBase);
931

Matthieu Schaller's avatar
Matthieu Schaller committed
932
933
        /* Number of updated particles */
        updated++;
934
        if (p->gpart != NULL) g_updated++;
935
936
937
      }

      /* Minimal time for next end of time-step */
938
939
      ti_end_min = min(p->ti_end, ti_end_min);
      ti_end_max = max(p->ti_end, ti_end_max);
940
    }
941
942
  }

943
  /* Otherwise, aggregate data from children. */
944
945
946
  else {

    /* Loop over the progeny. */
947
    for (int k = 0; k < 8; k++)