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

Pedro Gonnet's avatar
Pedro Gonnet committed
24
25
/* Config parameters. */
#include "../config.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
26
27
28
29

/* Some standard headers. */
#include <float.h>
#include <limits.h>
30
#include <stdlib.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
31

32
33
/* MPI headers. */
#ifdef WITH_MPI
34
#include <mpi.h>
35
36
#endif

37
38
39
/* This object's header. */
#include "runner.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
40
/* Local headers. */
41
#include "active.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
42
#include "approx_math.h"
43
#include "atomic.h"
44
#include "cell.h"
45
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
46
#include "cooling.h"
47
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
48
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
49
#include "engine.h"
50
#include "error.h"
51
52
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
53
#include "hydro_properties.h"
54
#include "kick.h"
55
#include "minmax.h"
56
#include "scheduler.h"
57
#include "sourceterms.h"
58
59
60
#include "space.h"
#include "task.h"
#include "timers.h"
61
#include "timestep.h"
62
#include "runner_doiact_vec.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
63

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

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

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

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

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

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

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

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

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

Tom Theuns's avatar
Tom Theuns committed
117
/**
Tom Theuns's avatar
Tom Theuns committed
118
 * @brief Perform source terms
Tom Theuns's avatar
Tom Theuns committed
119
120
121
122
123
124
125
 *
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
 */
void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
  const int count = c->count;
126
  const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
Tom Theuns's avatar
Tom Theuns committed
127
  const double cell_width[3] = {c->width[0], c->width[1], c->width[2]};
Tom Theuns's avatar
Tom Theuns committed
128
  struct sourceterms *sourceterms = r->e->sourceterms;
129
  const int dimen = 3;
Tom Theuns's avatar
Tom Theuns committed
130
131
132
133
134
135
136

  TIMER_TIC;

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_sourceterms(r, c->progeny[k], 0);
137
  } else {
Tom Theuns's avatar
Tom Theuns committed
138

139
    if (count > 0) {
Tom Theuns's avatar
Tom Theuns committed
140

141
142
143
144
145
146
      /* do sourceterms in this cell? */
      const int incell =
          sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
      if (incell == 1) {
        sourceterms_apply(r, sourceterms, c);
      }
Tom Theuns's avatar
Tom Theuns committed
147
148
    }
  }
Tom Theuns's avatar
Tom Theuns committed
149
150
151
152

  if (timer) TIMER_TOC(timer_dosource);
}

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

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

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

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

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

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

183
184
      /* Get a direct pointer on the part. */
      struct gpart *restrict gp = &gparts[i];
Matthieu Schaller's avatar
Matthieu Schaller committed
185

186
      /* Is this part within the time step? */
187
      if (gpart_is_active(gp, e)) {
188
189
        external_gravity_acceleration(time, potential, constants, gp);
      }
190
    }
191
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
192

193
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
194
195
}

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

  TIMER_TIC;

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
221
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
222

223
224
    /* Loop over the parts in this cell. */
    for (int i = 0; i < count; i++) {
Stefan Arridge's avatar
Stefan Arridge committed
225

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

230
231
      /* Kick has already updated ti_end, so need to check ti_begin */
      if (p->ti_begin == ti_current) {
232

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

235
236
        cooling_cool_part(constants, us, cooling_func, p, xp, dt);
      }
Stefan Arridge's avatar
Stefan Arridge committed
237
238
239
240
241
242
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
243
244
245
246
247
248
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
249
void runner_do_sort_ascending(struct entry *sort, int N) {
250
251
252
253
254
255
256
257
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

  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
304
        }
305
306
307
308
309
310
311
312
313
314
315
316
      } 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
317
    }
318
319
320
  }
}

Pedro Gonnet's avatar
Pedro Gonnet committed
321
322
323
324
325
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
326
 * @param flags Cell flag.
327
328
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
329
 */
330
void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
331
332
333
334
335
336
337

  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
338
339
  float buff[8];
  double px[3];
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363

  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;
364
      if (missing) runner_do_sort(r, c->progeny[k], missing, 0);
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
    }

    /* 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
417
        }
418

419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
      } /* 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
443
444
445
          sort[j * (count + 1) + k].d = px[0] * runner_shift[j][0] +
                                        px[1] * runner_shift[j][1] +
                                        px[2] * runner_shift[j][2];
446
        }
447
    }
448
449
450
451
452
453

    /* 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;
454
        runner_do_sort_ascending(&sort[j * (count + 1)], count);
455
456
457
458
        c->sorted |= (1 << j);
      }
  }

459
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
460
  /* Verify the sorting. */
461
462
463
464
465
466
467
468
469
470
  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
471
472
473
474

  if (clock) TIMER_TOC(timer_dosort);
}

475
476
477
478
479
/**
 * @brief Initialize the particles before the density calculation
 *
 * @param r The runner thread.
 * @param c The cell.
Matthieu Schaller's avatar
Matthieu Schaller committed
480
 * @param timer 1 if the time is to be recorded.
481
 */
482
void runner_do_init(struct runner *r, struct cell *c, int timer) {
483

Matthieu Schaller's avatar
Matthieu Schaller committed
484
485
  struct part *restrict parts = c->parts;
  struct gpart *restrict gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
486
  const int count = c->count;
487
  const int gcount = c->gcount;
488
  const struct engine *e = r->e;
489
490

  TIMER_TIC;
491

492
  /* Anything to do here? */
493
  if (!cell_is_active(c, e)) return;
494

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

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

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

507
      if (part_is_active(p, e)) {
Matthieu Schaller's avatar
Matthieu Schaller committed
508

509
510
        /* Get ready for a density calculation */
        hydro_init_part(p);
511
      }
512
    }
513
514
515
516
517

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

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

520
      if (gpart_is_active(gp, e)) {
521
522

        /* Get ready for a density calculation */
523
        gravity_init_gpart(gp);
524
525
      }
    }
526
  }
527

Peter W. Draper's avatar
Peter W. Draper committed
528
  if (timer) TIMER_TOC(timer_init);
529
530
}

531
/**
532
533
534
535
536
 * @brief Intermediate task after the gradient loop that does final operations
 * on the gradient quantities and optionally slope limits the gradients
 *
 * @param r The runner thread.
 * @param c The cell.
537
 * @param timer Are we timing this ?
538
 */
539
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
540

541
#ifdef EXTRA_HYDRO_LOOP
542

543
544
  struct part *restrict parts = c->parts;
  const int count = c->count;
545
  const struct engine *e = r->e;
546

547
548
  TIMER_TIC;

549
  /* Anything to do here? */
550
  if (!cell_is_active(c, e)) return;
551

552
553
554
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
555
      if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k], 0);
556
557
558
559
560
561
562
563
  } else {

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

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

564
      if (part_is_active(p, e)) {
565
566
567
568
569
570

        /* Get ready for a force calculation */
        hydro_end_gradient(p);
      }
    }
  }
571

572
573
  if (timer) TIMER_TOC(timer_do_extra_ghost);

574
575
#else
  error("SWIFT was not compiled with the extra hydro loop activated.");
576
#endif
577
}
578

579
/**
580
581
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
582
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
583
 * @param r The runner thread.
584
 * @param c The cell.
585
 * @param timer Are we timing this ?
586
 */
587
void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
588

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

602
603
  TIMER_TIC;

604
  /* Anything to do here? */
605
  if (!cell_is_active(c, e)) return;
606

607
608
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
609
    for (int k = 0; k < 8; k++)
610
611
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k], 0);
  } else {
612

613
614
615
616
617
    /* Init the IDs that have to be updated. */
    int *pid = NULL;
    if ((pid = malloc(sizeof(int) * count)) == NULL)
      error("Can't allocate memory for pid.");
    for (int k = 0; k < count; k++) pid[k] = k;
618

619
620
621
    /* While there are particles that need to be updated... */
    for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
         num_reruns++) {
622

623
624
      /* Reset the redo-count. */
      redo = 0;
625

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

629
630
631
        /* Get a direct pointer on the part. */
        struct part *restrict p = &parts[pid[i]];
        struct xpart *restrict xp = &xparts[pid[i]];
632

633
        /* Is this part within the timestep? */
634
        if (part_is_active(p, e)) {
635

636
          /* Finish the density calculation */
637
          hydro_end_density(p);
638

639
          float h_corr = 0.f;
Matthieu Schaller's avatar
Matthieu Schaller committed
640

641
642
          /* If no derivative, double the smoothing length. */
          if (p->density.wcount_dh == 0.0f) h_corr = p->h;
643

644
645
646
          /* Otherwise, compute the smoothing length update (Newton step). */
          else {
            h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
647

648
649
650
651
            /* Truncate to the range [ -p->h/2 , p->h ]. */
            h_corr = (h_corr < p->h) ? h_corr : p->h;
            h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
          }
652

653
654
655
          /* Did we get the right number density? */
          if (p->density.wcount > max_wcount ||
              p->density.wcount < min_wcount) {
656

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

660
661
662
            /* Flag for another round of fun */
            pid[redo] = pid[i];
            redo += 1;
663

664
665
            /* Re-initialise everything */
            hydro_init_part(p);
666

667
668
669
            /* Off we go ! */
            continue;
          }
670

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

673
          /* As of here, particle force variables will be set. */
674

675
676
          /* Compute variables required for the force loop */
          hydro_prepare_force(p, xp, ti_current, timeBase);
677

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

681
682
683
          /* Prepare the particle for the force loop over neighbours */
          hydro_reset_acceleration(p);
        }
684
685
      }

686
687
      /* We now need to treat the particles whose smoothing length had not
       * converged again */
688

689
690
691
      /* Re-set the counter for the next loop (potentially). */
      count = redo;
      if (count > 0) {
692

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

696
697
          /* 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
698

699
700
701
            /* Self-interaction? */
            if (l->t->type == task_type_self)
              runner_doself_subset_density(r, finger, parts, pid, count);
702

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

706
707
708
709
710
711
712
              /* 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);
713

714
            }
715

716
717
718
719
            /* 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);
720

721
722
723
724
725
726
727
728
729
730
731
            /* Otherwise, sub-pair interaction? */
            else if (l->t->type == task_type_sub_pair) {

              /* 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);
            }
732
733
734
          }
        }
      }
735
    }
736

737
738
    if (count)
      message("Smoothing length failed to converge on %i particles.", count);
739

740
741
742
    /* Be clean */
    free(pid);
  }
743

744
  if (timer) TIMER_TOC(timer_do_ghost);
745
746
}

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

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

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

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

781
782
783
  /* Now, we can drift */

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

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

795
796
  /* No children? */
  if (!c->split) {
797

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

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

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

808
809
        /* Drift... */
        drift_gpart(gp, dt, timeBase, ti_old, ti_current);
810

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

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

822
823
824
        /* Get a handle on the part. */
        struct part *const p = &parts[k];
        struct xpart *const xp = &xparts[k];
825

826
827
        /* Drift... */
        drift_part(p, xp, dt, timeBase, ti_old, ti_current);
Matthieu Schaller's avatar
Matthieu Schaller committed
828

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

835
836
837
        /* Maximal smoothing length */
        h_max = (h_max > p->h) ? h_max : p->h;
      }
838

839
840
      /* Now, get the maximal particle motion from its square */
      dx_max = sqrtf(dx2_max);
841

842
    } /* Check that we are actually going to move forward. */
843
844
845
846
847
848

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

Matthieu Schaller's avatar
Matthieu Schaller committed
851
852
853
  /* Otherwise, aggregate data from children. */
  else {

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

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

  /* Store the values */
867
868
869
870
871
  c->h_max = h_max;
  c->dx_max = dx_max;

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

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

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

885
886
887
888
889
  struct engine *e = (struct engine *)extra_data;
  struct cell *cells = (struct cell *)map_data;

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

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

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

917
  TIMER_TIC;
918

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

926
927
928
  int updated = 0, g_updated = 0;
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;

929
930
931
  /* No children? */
  if (!c->split) {

932
933
934
935
    /* Loop over the g-particles and kick the active ones. */
    for (int k = 0; k < gcount; k++) {

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

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

941
        if (gpart_is_active(gp, e)) {
942

Matthieu Schaller's avatar
Matthieu Schaller committed
943
          /* First, finish the force calculation */
Matthieu Schaller's avatar
Matthieu Schaller committed
944
          gravity_end_force(gp, const_G);
945

946
          /* Compute the next timestep */
947
          const int new_dti = get_gpart_timestep(gp, e);