runner.c 31.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 "const.h"
44
#include "debug.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
45
#include "engine.h"
46
#include "error.h"
47
#include "gravity.h"
48
#include "hydro_properties.h"
49
#include "hydro.h"
50
#include "drift.h"
51
#include "kick.h"
52
#include "minmax.h"
53
54
55
56
#include "scheduler.h"
#include "space.h"
#include "task.h"
#include "timers.h"
57
#include "timestep.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
58

59
/* Orientation of the cell pairs */
60
61
62
63
64
65
66
67
68
69
70
71
72
73
const float runner_shift[13 * 3] = {
    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,
    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, };
74
75

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

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

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

88
89
90
/* Import the gravity loop functions. */
#include "runner_doiact_grav.h"

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

100
101
102
  struct gpart *g, *gparts = c->gparts;
  int i, k, gcount = c->gcount;
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
103
  const struct external_potential *potential = r->e->external_potential;
104
  const struct phys_const *constants = r->e->physical_constants;
Matthieu Schaller's avatar
Matthieu Schaller committed
105

106
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
107
108
109
110

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

#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
116
  OUT;
Tom Theuns's avatar
Tom Theuns committed
117
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
118

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

122
123
124
125
126
    /* Get a direct pointer on the part. */
    g = &gparts[i];

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

128
      external_gravity(potential, constants, g);
129
    }
130
  }
131
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
132
133
}

Pedro Gonnet's avatar
Pedro Gonnet committed
134
135
136
137
138
139
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
140

141
void runner_do_sort_ascending(struct entry *sort, int N) {
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
185
186
187
188
189
190
191
192
193
194
195

  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
196
        }
197
198
199
200
201
202
203
204
205
206
207
208
      } 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
209
    }
210
211
212
  }
}

Pedro Gonnet's avatar
Pedro Gonnet committed
213
214
215
216
217
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
218
 * @param flags Cell flag.
219
220
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
221
 */
222

223
void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256

  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;
  // float shift[3];
  float buff[8], px[3];

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

    /* 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
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
      } /* 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;
          sort[j * (count + 1) + k].d = px[0] * runner_shift[3 * j + 0] +
                                        px[1] * runner_shift[3 * j + 1] +
                                        px[2] * runner_shift[3 * j + 2];
        }
340
    }
341
342
343
344
345
346

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

Matthieu Schaller's avatar
Matthieu Schaller committed
352
353
354
355
356
357
358
359
360
361
362
363
  /* Verify the sorting. */
  /* 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." );
          }
      } */
364
365
366
367

  if (clock) TIMER_TOC(timer_dosort);
}

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

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

  TIMER_TIC;
384

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

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

395
      /* Get a direct pointer on the part. */
396
      struct part *const p = &parts[i];
397

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

400
401
        /* Get ready for a density calculation */
        hydro_init_part(p);
402
      }
403
    }
404
405
406
407
408
409
410
411
412
413
414

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

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

      if (gp->ti_end <= ti_current) {

        /* Get ready for a density calculation */
        gravity_init_part(gp);
415
416
417
418

        if (gp->id == -ICHECK)
          message("id=%lld a=[%f %f %f]\n", gp->id, gp->a_grav[0],
                  gp->a_grav[1], gp->a_grav[2]);
419
420
      }
    }
421
  }
422

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

426
427
428
/**
 * @brief Intermediate task between density and force
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
429
 * @param r The runner thread.
430
 * @param c The cell.
431
 */
432

Matthieu Schaller's avatar
Matthieu Schaller committed
433
void runner_do_ghost(struct runner *r, struct cell *c) {
434
435

  struct part *p, *parts = c->parts;
436
  struct xpart *xp, *xparts = c->xparts;
437
  struct cell *finger;
Matthieu Schaller's avatar
Matthieu Schaller committed
438
  int redo, count = c->count;
439
  int *pid;
440
  float h_corr;
441
442
  const int ti_current = r->e->ti_current;
  const double timeBase = r->e->timeBase;
443
444
445
446
447
448
449
  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;
450

451
452
  TIMER_TIC;

453
454
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
455
    for (int k = 0; k < 8; k++)
Matthieu Schaller's avatar
Matthieu Schaller committed
456
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k]);
457
458
459
460
461
462
    return;
  }

  /* Init the IDs that have to be updated. */
  if ((pid = (int *)alloca(sizeof(int) * count)) == NULL)
    error("Call to alloca failed.");
Matthieu Schaller's avatar
Matthieu Schaller committed
463
  for (int k = 0; k < count; k++) pid[k] = k;
464
465

  /* While there are particles that need to be updated... */
466
  for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
Matthieu Schaller's avatar
Matthieu Schaller committed
467
       num_reruns++) {
468
469
470

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

472
    /* Loop over the parts in this cell. */
Matthieu Schaller's avatar
Matthieu Schaller committed
473
    for (int i = 0; i < count; i++) {
474
475
476

      /* Get a direct pointer on the part. */
      p = &parts[pid[i]];
477
      xp = &xparts[pid[i]];
478
479

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

482
        /* Finish the density calculation */
483
        hydro_end_density(p, ti_current);
484
485

        /* If no derivative, double the smoothing length. */
486
        if (p->density.wcount_dh == 0.0f) h_corr = p->h;
487
488
489

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

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

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

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

503
          /* Flag for another round of fun */
504
505
          pid[redo] = pid[i];
          redo += 1;
506

507
508
          /* Re-initialise everything */
          hydro_init_part(p);
509

510
          /* Off we go ! */
511
          continue;
512
513
        }

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

516
        /* As of here, particle force variables will be set. */
517

518
        /* Compute variables required for the force loop */
519
        hydro_prepare_force(p, xp, ti_current, timeBase);
520

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

524
525
        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
526
527
528
      }
    }

529
530
531
    /* We now need to treat the particles whose smoothing length had not
     * converged again */

532
533
534
535
536
537
    /* Re-set the counter for the next loop (potentially). */
    count = redo;
    if (count > 0) {

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

539
540
        /* 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
541

542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
          /* 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);

          }

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

            /* 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);
          }
        }
      }
572
    }
573
  }
574

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

Matthieu Schaller's avatar
Matthieu Schaller committed
578
  TIMER_TOC(timer_do_ghost);
579
580
}

581
/**
582
 * @brief Drift particles and g-particles forward in time
583
584
585
 *
 * @param r The runner thread.
 * @param c The cell.
586
 * @param timer Are we timing this ?
587
 */
588
void runner_do_drift(struct runner *r, struct cell *c, int timer) {
589

590
591
  const double timeBase = r->e->timeBase;
  const double dt = (r->e->ti_current - r->e->ti_old) * timeBase;
592
593
594
595
596
  const int ti_old = r->e->ti_old;
  const int ti_current = r->e->ti_current;
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;
597
  float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
598

599
600
601
602
  double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
  double mom[3] = {0.0, 0.0, 0.0};
  double ang_mom[3] = {0.0, 0.0, 0.0};

603
  TIMER_TIC
604

Tom Theuns's avatar
Tom Theuns committed
605
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
606
  OUT;
Tom Theuns's avatar
Tom Theuns committed
607
#endif
608

609
610
  /* No children? */
  if (!c->split) {
611

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

616
617
618
619
      /* Get a handle on the gpart. */
      struct gpart *const gp = &gparts[k];

      /* Drift... */
620
      drift_gpart(gp, dt, timeBase, ti_old, ti_current);
621
622
623
624
625
626

      /* 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);
627
628
629
    }

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

633
      /* Get a handle on the part. */
634
635
      struct part *const p = &parts[k];
      struct xpart *const xp = &xparts[k];
636

Matthieu Schaller's avatar
Matthieu Schaller committed
637

638
      /* Drift... */
639
      drift_part(p, xp, dt, timeBase, ti_old, ti_current);
640

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

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

      /* 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.;
676
      e_int += m * hydro_get_internal_energy(p, half_dt);
677
    }
678
679

    /* Now, get the maximal particle motion from its square */
680
    dx_max = sqrtf(dx2_max);
681
  }
682

Matthieu Schaller's avatar
Matthieu Schaller committed
683
684
685
686
687
688
  /* Otherwise, aggregate data from children. */
  else {

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

        /* Recurse */
Matthieu Schaller's avatar
Matthieu Schaller committed
691
        struct cell *cp = c->progeny[k];
692
        runner_do_drift(r, cp, 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
693

694
        /* Collect */
695
696
        dx_max = fmaxf(dx_max, cp->dx_max);
        h_max = fmaxf(h_max, cp->h_max);
697
698
699
700
701
702
703
704
705
706
        mass += cp->mass;
        e_kin += cp->e_kin;
        e_int += cp->e_int;
        e_pot += cp->e_pot;
        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
707
708
709
710
711
712
      }
  }

  /* Store the values */
  c->h_max = h_max;
  c->dx_max = dx_max;
713
714
715
716
717
718
719
720
721
722
  c->mass = mass;
  c->e_kin = e_kin;
  c->e_int = e_int;
  c->e_pot = e_pot;
  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];
723

Peter W. Draper's avatar
Peter W. Draper committed
724
  if (timer) TIMER_TOC(timer_drift);
725
}
726

727
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
728
729
 * @brief Kick particles in momentum space and collect statistics (fixed
 * time-step case)
730
731
732
 *
 * @param r The runner thread.
 * @param c The cell.
733
 * @param timer Are we timing this ?
734
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
735
void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
736

Matthieu Schaller's avatar
Matthieu Schaller committed
737
  const double global_dt = r->e->dt_max;
738
  const double timeBase = r->e->timeBase;
Matthieu Schaller's avatar
Matthieu Schaller committed
739
  const int count = c->count;
740
741
742
743
  const int gcount = c->gcount;
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
744
  const double const_G = r->e->physical_constants->const_newton_G;
Matthieu Schaller's avatar
Matthieu Schaller committed
745

746
  int updated = 0, g_updated = 0;
747
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
748
749
750

  TIMER_TIC

Tom Theuns's avatar
Tom Theuns committed
751
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
752
  OUT;
Tom Theuns's avatar
Tom Theuns committed
753
754
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
755
756
757
  /* The new time-step */
  const int new_dti = global_dt / timeBase;

758
759
760
  /* No children? */
  if (!c->split) {

761
    /* Loop over the g-particles and kick everyone. */
762
763
764
765
766
    for (int k = 0; k < gcount; k++) {

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

767
      /* If the g-particle has no counterpart */
768
      if (gp->id < 0) {
769

Matthieu Schaller's avatar
Matthieu Schaller committed
770
        /* First, finish the force calculation */
771
        gravity_end_force(gp, const_G);
772

773
774
        /* Kick the g-particle forward */
        kick_gpart(gp, new_dti, timeBase);
775

Matthieu Schaller's avatar
Matthieu Schaller committed
776
777
        /* Number of updated g-particles */
        g_updated++;
778

Matthieu Schaller's avatar
Matthieu Schaller committed
779
780
781
782
783
        /* 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);
      }
    }
784

Matthieu Schaller's avatar
Matthieu Schaller committed
785
    /* Now do the hydro ones... */
786

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

Matthieu Schaller's avatar
Matthieu Schaller committed
790
791
792
      /* Get a handle on the part. */
      struct part *const p = &parts[k];
      struct xpart *const xp = &xparts[k];
793

Matthieu Schaller's avatar
Matthieu Schaller committed
794
795
      /* First, finish the force loop */
      p->h_dt *= p->h * 0.333333333f;
796

Matthieu Schaller's avatar
Matthieu Schaller committed
797
798
      /* And do the same of the extra variable */
      hydro_end_force(p);
799
      if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
800

801
802
      /* Kick the particle forward */
      kick_part(p, xp, new_dti, timeBase);
803

Matthieu Schaller's avatar
Matthieu Schaller committed
804
805
806
      /* Number of updated particles */
      updated++;
      if (p->gpart != NULL) g_updated++;
807

Matthieu Schaller's avatar
Matthieu Schaller committed
808
809
810
      /* 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
811
812
    }
  }
813

Matthieu Schaller's avatar
Matthieu Schaller committed
814
815
  /* Otherwise, aggregate data from children. */
  else {
816

Matthieu Schaller's avatar
Matthieu Schaller committed
817
818
819
820
    /* Loop over the progeny. */
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
        struct cell *const cp = c->progeny[k];
821

Matthieu Schaller's avatar
Matthieu Schaller committed
822
823
        /* Recurse */
        runner_do_kick_fixdt(r, cp, 0);
824

Matthieu Schaller's avatar
Matthieu Schaller committed
825
826
827
828
829
830
831
        /* 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);
      }
  }
832

Matthieu Schaller's avatar
Matthieu Schaller committed
833
834
835
836
837
  /* 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;
838

Matthieu Schaller's avatar
Matthieu Schaller committed
839
840
  if (timer) TIMER_TOC(timer_kick);
}
841

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

852
853
  const struct engine *e = r->e;
  const double timeBase = e->timeBase;
854
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
855
  const int count = c->count;
856
857
858
859
  const int gcount = c->gcount;
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
860
  const double const_G = r->e->physical_constants->const_newton_G;
Matthieu Schaller's avatar
Matthieu Schaller committed
861

862
  int updated = 0, g_updated = 0;
863
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
864
865

  TIMER_TIC
866

Tom Theuns's avatar
Tom Theuns committed
867
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
868
  OUT;
Tom Theuns's avatar
Tom Theuns committed
869
870
#endif

871
872
873
  /* No children? */
  if (!c->split) {

874
875
876
877
878
879
880
    /* Loop over the g-particles and kick the active ones. */
    for (int k = 0; k < gcount; k++) {

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

      /* If the g-particle has no counterpart and needs to be kicked */
881
      if (gp->id < 0) {
882

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

Matthieu Schaller's avatar
Matthieu Schaller committed
885
          /* First, finish the force calculation */
Matthieu Schaller's avatar
Matthieu Schaller committed
886
          gravity_end_force(gp, const_G);
887

888
          /* Compute the next timestep */
889
          const int new_dti = get_gpart_timestep(gp, e);
890

Matthieu Schaller's avatar
Matthieu Schaller committed
891
          /* Now we have a time step, proceed with the kick */
892
          kick_gpart(gp, new_dti, timeBase);
Matthieu Schaller's avatar
Matthieu Schaller committed
893
894
895
896
897
898
899
900

          /* 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);
901
902
903
904
905
      }
    }

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

906
    /* Loop over the particles and kick the active ones. */
Matthieu Schaller's avatar
Matthieu Schaller committed
907
    for (int k = 0; k < count; k++) {
908
909

      /* Get a handle on the part. */
Matthieu Schaller's avatar
Matthieu Schaller committed
910
911
      struct part *const p = &parts[k];
      struct xpart *const xp = &xparts[k];
912
913

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

        /* First, finish the force loop */
917
        p->h_dt *= p->h * 0.333333333f;
918
919

        /* And do the same of the extra variable */
920
        hydro_end_force(p);
Matthieu Schaller's avatar
Matthieu Schaller committed
921
        if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
922

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

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

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

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

940
  /* Otherwise, aggregate data from children. */
941
942
943
  else {

    /* Loop over the progeny. */
944
    for (int k = 0; k < 8; k++)
945
      if (c->progeny[k] != NULL) {
Matthieu Schaller's avatar
Matthieu Schaller committed
946
        struct cell *const cp = c->progeny[k];
947
948

        /* Recurse */
949
        runner_do_kick(r, cp, 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
950

951
        /* And aggregate */
952
        updated += cp->updated;
953
        g_updated += cp->g_updated;
954
955
        ti_end_min = min(cp->ti_end_min, ti_end_min);
        ti_end_max = max(cp->ti_end_max, ti_end_max);
956
957
958
959
      }
  }

  /* Store the values. */
960
  c->updated = updated;
961
  c->g_updated = g_updated;
962
963
  c->ti_end_min = ti_end_min;
  c->ti_end_max = ti_end_max;
964

Peter W. Draper's avatar
Peter W. Draper committed
965
  if (timer) TIMER_TOC(timer_kick);
966
}
967

968
969
970
971
972
973
974
/**
 * @brief Construct the cell properties from the received particles
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
975
void runner_do_recv_cell(struct runner *r, struct cell *c, int timer) {
976
977
978
979
980

  const struct part *const parts = c->parts;
  const struct gpart *const gparts = c->gparts;
  const size_t nr_parts = c->count;
  const size_t nr_gparts = c->gcount;
Matthieu Schaller's avatar
Matthieu Schaller committed
981
  // const int ti_current = r->e->ti_current;
982
983
984
985
986
987
988

  TIMER_TIC;

  int ti_end_min = max_nr_timesteps;
  int ti_end_max = 0;
  float h_max = 0.f;

989
  /* Collect everything... */
990