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

Stefan Arridge's avatar
Stefan Arridge committed
138
139
140
141
142
143
144
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

/**
 * @brief Calculate change in entropy from cooling
 *
 * @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;
  const int count = c->count;
  const int ti_current = r->e->ti_current;
  const struct cooling_data *cooling = r->e->cooling;
  const struct phys_const *constants = r->e->physical_constants;
  const double timeBase = r->e->timeBase;
  double dt;

  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

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

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

    /* Is this part within the time step? */
    if (p->ti_end <= ti_current) {
      dt = (p->ti_end - p->ti_begin)*timeBase;
      update_entropy(cooling, constants, p, dt);
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
185
186
187
188
189
190
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
191
void runner_do_sort_ascending(struct entry *sort, int N) {
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245

  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
246
        }
247
248
249
250
251
252
253
254
255
256
257
258
      } 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
259
    }
260
261
262
  }
}

Pedro Gonnet's avatar
Pedro Gonnet committed
263
264
265
266
267
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
268
 * @param flags Cell flag.
269
270
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
271
 */
272
void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
273
274
275
276
277
278
279

  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
280
281
  float buff[8];
  double px[3];
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305

  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;
306
      if (missing) runner_do_sort(r, c->progeny[k], missing, 0);
307
308
309
310
311
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
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
    }

    /* 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
359
        }
360

361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
      } /* 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
385
386
387
          sort[j * (count + 1) + k].d = px[0] * runner_shift[j][0] +
                                        px[1] * runner_shift[j][1] +
                                        px[2] * runner_shift[j][2];
388
        }
389
    }
390
391
392
393
394
395

    /* 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;
396
        runner_do_sort_ascending(&sort[j * (count + 1)], count);
397
398
399
400
        c->sorted |= (1 << j);
      }
  }

401
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
402
  /* Verify the sorting. */
403
404
405
406
407
408
409
410
411
412
  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
413
414
415
416

  if (clock) TIMER_TOC(timer_dosort);
}

417
418
419
420
421
/**
 * @brief Initialize the particles before the density calculation
 *
 * @param r The runner thread.
 * @param c The cell.
Matthieu Schaller's avatar
Matthieu Schaller committed
422
 * @param timer 1 if the time is to be recorded.
423
 */
424
void runner_do_init(struct runner *r, struct cell *c, int timer) {
425

Matthieu Schaller's avatar
Matthieu Schaller committed
426
427
  struct part *restrict parts = c->parts;
  struct gpart *restrict gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
428
  const int count = c->count;
429
  const int gcount = c->gcount;
430
  const int ti_current = r->e->ti_current;
431
432

  TIMER_TIC;
433

434
435
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
436
    for (int k = 0; k < 8; k++)
437
      if (c->progeny[k] != NULL) runner_do_init(r, c->progeny[k], 0);
438
    return;
439
440
  } else {

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

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

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

449
450
        /* Get ready for a density calculation */
        hydro_init_part(p);
451
      }
452
    }
453
454
455
456
457

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

      /* Get a direct pointer on the part. */
458
      struct gpart *restrict gp = &gparts[i];
459
460
461
462

      if (gp->ti_end <= ti_current) {

        /* Get ready for a density calculation */
463
        gravity_init_gpart(gp);
464
465
      }
    }
466
  }
467

Peter W. Draper's avatar
Peter W. Draper committed
468
  if (timer) TIMER_TOC(timer_init);
469
470
}

471
472
473
/**
 * @brief Intermediate task between density and force
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
474
 * @param r The runner thread.
475
 * @param c The cell.
476
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
477
void runner_do_ghost(struct runner *r, struct cell *c) {
478

Matthieu Schaller's avatar
Matthieu Schaller committed
479
480
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
481
  int redo, count = c->count;
482
483
  const int ti_current = r->e->ti_current;
  const double timeBase = r->e->timeBase;
484
485
486
487
488
489
490
  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;
491

492
493
  TIMER_TIC;

494
495
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
496
    for (int k = 0; k < 8; k++)
Matthieu Schaller's avatar
Matthieu Schaller committed
497
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k]);
498
499
500
501
    return;
  }

  /* Init the IDs that have to be updated. */
502
503
504
  int *pid = NULL;
  if ((pid = malloc(sizeof(int) * count)) == NULL)
    error("Can't allocate memory for pid.");
Matthieu Schaller's avatar
Matthieu Schaller committed
505
  for (int k = 0; k < count; k++) pid[k] = k;
506
507

  /* While there are particles that need to be updated... */
508
  for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
Matthieu Schaller's avatar
Matthieu Schaller committed
509
       num_reruns++) {
510
511
512

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

514
    /* Loop over the parts in this cell. */
Matthieu Schaller's avatar
Matthieu Schaller committed
515
    for (int i = 0; i < count; i++) {
516
517

      /* Get a direct pointer on the part. */
518
519
      struct part *restrict p = &parts[pid[i]];
      struct xpart *restrict xp = &xparts[pid[i]];
520
521

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

524
        /* Finish the density calculation */
525
        hydro_end_density(p, ti_current);
526

Matthieu Schaller's avatar
Matthieu Schaller committed
527
528
        float h_corr = 0.f;

529
        /* If no derivative, double the smoothing length. */
530
        if (p->density.wcount_dh == 0.0f) h_corr = p->h;
531
532
533

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

          /* Truncate to the range [ -p->h/2 , p->h ]. */
537
538
          h_corr = fminf(h_corr, p->h);
          h_corr = fmaxf(h_corr, -p->h * 0.5f);
Pedro Gonnet's avatar
Pedro Gonnet committed
539
        }
540
541

        /* Did we get the right number density? */
542
        if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
543
544
545
546

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

547
          /* Flag for another round of fun */
548
549
          pid[redo] = pid[i];
          redo += 1;
550

551
552
          /* Re-initialise everything */
          hydro_init_part(p);
553

554
          /* Off we go ! */
555
          continue;
556
557
        }

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

560
        /* As of here, particle force variables will be set. */
561

562
        /* Compute variables required for the force loop */
563
        hydro_prepare_force(p, xp, ti_current, timeBase);
564

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

568
569
        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
570
571
572
      }
    }

573
574
575
    /* We now need to treat the particles whose smoothing length had not
     * converged again */

576
577
578
579
580
    /* 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
581
      for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
Matthieu Schaller's avatar
Matthieu Schaller committed
582

583
584
        /* 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
585

586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
          /* 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);

          }

603
604
605
606
607
608
609
          /* 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) {
610
611
612
613
614
615
616
617
618
619
620

            /* 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);
          }
        }
      }
621
    }
622
  }
623

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

627
628
629
  /* Be clean */
  free(pid);

Matthieu Schaller's avatar
Matthieu Schaller committed
630
  TIMER_TOC(timer_do_ghost);
631
632
}

633
/**
634
 * @brief Drift particles and g-particles forward in time
635
636
637
 *
 * @param r The runner thread.
 * @param c The cell.
638
 * @param timer Are we timing this ?
639
 */
640
void runner_do_drift(struct runner *r, struct cell *c, int timer) {
641

642
643
  const double timeBase = r->e->timeBase;
  const double dt = (r->e->ti_current - r->e->ti_old) * timeBase;
644
645
  const int ti_old = r->e->ti_old;
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
646
647
648
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
  struct gpart *restrict gparts = c->gparts;
649
  float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
650

651
  double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
652
653
654
  double mom[3] = {0.0, 0.0, 0.0};
  double ang_mom[3] = {0.0, 0.0, 0.0};

655
  TIMER_TIC
656

Tom Theuns's avatar
Tom Theuns committed
657
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
658
  OUT;
Tom Theuns's avatar
Tom Theuns committed
659
#endif
660

661
662
  /* No children? */
  if (!c->split) {
663

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

668
      /* Get a handle on the gpart. */
669
      struct gpart *restrict gp = &gparts[k];
670
671

      /* Drift... */
672
      drift_gpart(gp, dt, timeBase, ti_old, ti_current);
673
674
675
676
677
678

      /* 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);
679
680
681
    }

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

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

689
      /* Drift... */
690
      drift_part(p, xp, dt, timeBase, ti_old, ti_current);
691

692
      /* Compute (square of) motion since last cell construction */
693
694
695
      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];
696
      dx2_max = fmaxf(dx2_max, dx2);
Matthieu Schaller's avatar
Matthieu Schaller committed
697
698
699

      /* Maximal smoothing length */
      h_max = fmaxf(p->h, h_max);
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726

      /* 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.;
727
      e_int += m * hydro_get_internal_energy(p, half_dt);
728
729

      /* Collect entropy */
730
      entropy += m * hydro_get_entropy(p, half_dt);
731
    }
732
733

    /* Now, get the maximal particle motion from its square */
734
    dx_max = sqrtf(dx2_max);
735
  }
736

Matthieu Schaller's avatar
Matthieu Schaller committed
737
738
739
740
741
742
  /* Otherwise, aggregate data from children. */
  else {

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

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

748
        /* Collect */
749
750
        dx_max = fmaxf(dx_max, cp->dx_max);
        h_max = fmaxf(h_max, cp->h_max);
751
752
753
754
        mass += cp->mass;
        e_kin += cp->e_kin;
        e_int += cp->e_int;
        e_pot += cp->e_pot;
755
        entropy += cp->entropy;
756
757
758
759
760
761
        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
762
763
764
765
766
767
      }
  }

  /* Store the values */
  c->h_max = h_max;
  c->dx_max = dx_max;
768
769
770
771
  c->mass = mass;
  c->e_kin = e_kin;
  c->e_int = e_int;
  c->e_pot = e_pot;
772
  c->entropy = entropy;
773
774
775
776
777
778
  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];
779

Peter W. Draper's avatar
Peter W. Draper committed
780
  if (timer) TIMER_TOC(timer_drift);
781
}
782

783
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
784
785
 * @brief Kick particles in momentum space and collect statistics (fixed
 * time-step case)
786
787
788
 *
 * @param r The runner thread.
 * @param c The cell.
789
 * @param timer Are we timing this ?
790
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
791
void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
792

Matthieu Schaller's avatar
Matthieu Schaller committed
793
  const double global_dt = r->e->dt_max;
794
  const double timeBase = r->e->timeBase;
Matthieu Schaller's avatar
Matthieu Schaller committed
795
  const int count = c->count;
796
  const int gcount = c->gcount;
Matthieu Schaller's avatar
Matthieu Schaller committed
797
798
799
  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
800
  const double const_G = r->e->physical_constants->const_newton_G;
Matthieu Schaller's avatar
Matthieu Schaller committed
801

802
  int updated = 0, g_updated = 0;
803
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
804
805
806

  TIMER_TIC

Tom Theuns's avatar
Tom Theuns committed
807
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
808
  OUT;
Tom Theuns's avatar
Tom Theuns committed
809
810
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
811
812
813
  /* The new time-step */
  const int new_dti = global_dt / timeBase;

814
815
816
  /* No children? */
  if (!c->split) {

817
    /* Loop over the g-particles and kick everyone. */
818
819
820
    for (int k = 0; k < gcount; k++) {

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

823
      /* If the g-particle has no counterpart */
824
      if (gp->id_or_neg_offset > 0) {
825

Matthieu Schaller's avatar
Matthieu Schaller committed
826
        /* First, finish the force calculation */
827
        gravity_end_force(gp, const_G);
828

829
830
        /* Kick the g-particle forward */
        kick_gpart(gp, new_dti, timeBase);
831

Matthieu Schaller's avatar
Matthieu Schaller committed
832
833
        /* Number of updated g-particles */
        g_updated++;
834

Matthieu Schaller's avatar
Matthieu Schaller committed
835
836
837
838
839
        /* 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);
      }
    }
840

Matthieu Schaller's avatar
Matthieu Schaller committed
841
    /* Now do the hydro ones... */
842

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

Matthieu Schaller's avatar
Matthieu Schaller committed
846
      /* Get a handle on the part. */
847
848
      struct part *restrict p = &parts[k];
      struct xpart *restrict xp = &xparts[k];
849

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

854
855
      /* Kick the particle forward */
      kick_part(p, xp, new_dti, timeBase);
856

Matthieu Schaller's avatar
Matthieu Schaller committed
857
858
859
      /* Number of updated particles */
      updated++;
      if (p->gpart != NULL) g_updated++;
860

Matthieu Schaller's avatar
Matthieu Schaller committed
861
862
863
      /* 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
864
865
    }
  }
866

Matthieu Schaller's avatar
Matthieu Schaller committed
867
868
  /* Otherwise, aggregate data from children. */
  else {
869

Matthieu Schaller's avatar
Matthieu Schaller committed
870
871
872
    /* Loop over the progeny. */
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
873
        struct cell *restrict cp = c->progeny[k];
874

Matthieu Schaller's avatar
Matthieu Schaller committed
875
876
        /* Recurse */
        runner_do_kick_fixdt(r, cp, 0);
877

Matthieu Schaller's avatar
Matthieu Schaller committed
878
879
880
881
882
883
884
        /* 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);
      }
  }
885

Matthieu Schaller's avatar
Matthieu Schaller committed
886
887
888
889
890
  /* 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;
891

Matthieu Schaller's avatar
Matthieu Schaller committed
892
893
  if (timer) TIMER_TOC(timer_kick);
}
894

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

905
906
  const struct engine *e = r->e;
  const double timeBase = e->timeBase;
907
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
908
  const int count = c->count;
909
  const int gcount = c->gcount;
Matthieu Schaller's avatar
Matthieu Schaller committed
910
911
912
  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
913
  const double const_G = r->e->physical_constants->const_newton_G;
Matthieu Schaller's avatar
Matthieu Schaller committed
914

915
  int updated = 0, g_updated = 0;
916
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
917
918

  TIMER_TIC
919

Tom Theuns's avatar
Tom Theuns committed
920
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
921
  OUT;
Tom Theuns's avatar
Tom Theuns committed
922
923
#endif

924
925
926
  /* No children? */
  if (!c->split) {

927
928
929
930
    /* Loop over the g-particles and kick the active ones. */
    for (int k = 0; k < gcount; k++) {

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

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
938
          /* First, finish the force calculation */
Matthieu Schaller's avatar
Matthieu Schaller committed
939
          gravity_end_force(gp, const_G);
940

941
          /* Compute the next timestep */
942
          const int new_dti = get_gpart_timestep(gp, e);
943

Matthieu Schaller's avatar
Matthieu Schaller committed
944
          /* Now we have a time step, proceed with the kick */
945
          kick_gpart(gp, new_dti, timeBase);
Matthieu Schaller's avatar
Matthieu Schaller committed
946
947
948
949
950
951
952
953

          /* 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);
954
955
956
957
958
      }
    }

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

959
    /* Loop over the particles and kick the active ones. */
Matthieu Schaller's avatar
Matthieu Schaller committed
960
    for (int k = 0; k < count; k++) {
961
962

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

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

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

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

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