runner.c 49.3 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"
James Willis's avatar
James Willis committed
56
#include "runner_doiact_vec.h"
57
#include "scheduler.h"
58
#include "sourceterms.h"
59
#include "space.h"
60
#include "stars.h"
61
62
#include "task.h"
#include "timers.h"
63
#include "timestep.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
64

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

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

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

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

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

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

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

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

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

Tom Theuns's avatar
Tom Theuns committed
118
/**
Tom Theuns's avatar
Tom Theuns committed
119
 * @brief Perform source terms
Tom Theuns's avatar
Tom Theuns committed
120
121
122
123
124
125
126
 *
 * @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;
127
  const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
Tom Theuns's avatar
Tom Theuns committed
128
  const double cell_width[3] = {c->width[0], c->width[1], c->width[2]};
Tom Theuns's avatar
Tom Theuns committed
129
  struct sourceterms *sourceterms = r->e->sourceterms;
130
  const int dimen = 3;
Tom Theuns's avatar
Tom Theuns committed
131
132
133
134
135
136
137

  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);
138
  } else {
Tom Theuns's avatar
Tom Theuns committed
139

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

142
143
144
145
146
147
      /* 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
148
149
    }
  }
Tom Theuns's avatar
Tom Theuns committed
150
151
152
153

  if (timer) TIMER_TOC(timer_dosource);
}

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

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

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

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

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

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

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

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

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

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

  TIMER_TIC;

218
219
220
  /* Anything to do here? */
  if (!cell_is_active(c, e)) return;

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

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

230
231
232
      /* 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
233

234
      if (part_is_active(p, e)) {
235

236
237
        /* Let's cool ! */
        const double dt = get_timestep(p->time_bin, timeBase);
238
239
        cooling_cool_part(constants, us, cooling_func, p, xp, dt);
      }
Stefan Arridge's avatar
Stefan Arridge committed
240
241
242
243
244
245
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

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

  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
307
        }
308
309
310
311
312
313
314
315
316
317
318
319
      } 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
320
    }
321
322
323
  }
}

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

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

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

    /* 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
420
        }
421

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

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

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

  if (clock) TIMER_TOC(timer_dosort);
}

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

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

  TIMER_TIC;
494

495
  /* Anything to do here? */
496
  if (!cell_is_active(c, e)) return;
497

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

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

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

510
      if (part_is_active(p, e)) {
Matthieu Schaller's avatar
Matthieu Schaller committed
511

512
513
        /* Get ready for a density calculation */
        hydro_init_part(p);
514
      }
515
    }
516
517
518
519
520

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

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

523
      if (gpart_is_active(gp, e)) {
524
525

        /* Get ready for a density calculation */
526
        gravity_init_gpart(gp);
527
528
      }
    }
529
  }
530

Peter W. Draper's avatar
Peter W. Draper committed
531
  if (timer) TIMER_TOC(timer_init);
532
533
}

534
/**
535
536
537
538
539
 * @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.
540
 * @param timer Are we timing this ?
541
 */
542
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
543

544
#ifdef EXTRA_HYDRO_LOOP
545

546
547
  struct part *restrict parts = c->parts;
  const int count = c->count;
548
  const struct engine *e = r->e;
549

550
551
  TIMER_TIC;

552
  /* Anything to do here? */
553
  if (!cell_is_active(c, e)) return;
554

555
556
557
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
558
      if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k], 0);
559
560
561
562
563
564
565
566
  } 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];

567
      if (part_is_active(p, e)) {
568
569
570
571
572
573

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

575
576
  if (timer) TIMER_TOC(timer_do_extra_ghost);

577
578
#else
  error("SWIFT was not compiled with the extra hydro loop activated.");
579
#endif
580
}
581

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

Matthieu Schaller's avatar
Matthieu Schaller committed
592
593
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
594
595
  const struct engine *e = r->e;
  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
  int redo = 0, count = 0;
602

603
604
  TIMER_TIC;

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

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

614
    /* Init the list of active particles that have to be updated. */
615
    int *pid = NULL;
616
    if ((pid = malloc(sizeof(int) * c->count)) == NULL)
617
      error("Can't allocate memory for pid.");
618
619
620
621
622
    for (int k = 0; k < c->count; k++)
      if (part_is_active(&parts[k], e)) {
        pid[count] = k;
        ++count;
      }
623

624
625
626
    /* While there are particles that need to be updated... */
    for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
         num_reruns++) {
627

628
629
      /* Reset the redo-count. */
      redo = 0;
630

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

634
635
636
        /* Get a direct pointer on the part. */
        struct part *restrict p = &parts[pid[i]];
        struct xpart *restrict xp = &xparts[pid[i]];
637

638
#ifdef SWIFT_DEBUG_CHECKS
639
        /* Is this part within the timestep? */
640
641
642
643
644
        if (!part_is_active(p, e)) error("Ghost applied to inactive particle");
#endif

        /* Finish the density calculation */
        hydro_end_density(p);
645

646
647
        /* Did we get the right number of neighbours? */
        if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
648

649
          float h_corr = 0.f;
Matthieu Schaller's avatar
Matthieu Schaller committed
650

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

654
655
656
          /* Otherwise, compute the smoothing length update (Newton step). */
          else {
            h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
657

658
659
660
661
            /* 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;
          }
662

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

666
667
668
          /* Flag for another round of fun */
          pid[redo] = pid[i];
          redo += 1;
669

670
671
          /* Re-initialise everything */
          hydro_init_part(p);
672

673
674
675
          /* Off we go ! */
          continue;
        }
676

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

679
        /* As of here, particle force variables will be set. */
680

681
682
        /* Compute variables required for the force loop */
        hydro_prepare_force(p, xp);
683

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

687
688
        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
689
690
      }

691
692
      /* We now need to treat the particles whose smoothing length had not
       * converged again */
693

694
695
696
      /* Re-set the counter for the next loop (potentially). */
      count = redo;
      if (count > 0) {
697

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

701
702
          /* 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
703

704
705
706
            /* Self-interaction? */
            if (l->t->type == task_type_self)
              runner_doself_subset_density(r, finger, parts, pid, count);
707

708
709
            /* Otherwise, pair interaction? */
            else if (l->t->type == task_type_pair) {
710

711
712
713
714
715
716
717
              /* 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);
718

719
            }
720

721
722
723
724
            /* 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);
725

726
727
728
729
730
731
732
733
734
735
736
            /* 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);
            }
737
738
739
          }
        }
      }
740
    }
741

742
743
744
745
746
747
748
#ifdef SWIFT_DEBUG_CHECKS
    if (count) {
      message("Smoothing length failed to converge on %i particles.", count);

      error("Aborting....");
    }
#else
749
750
    if (count)
      message("Smoothing length failed to converge on %i particles.", count);
751
#endif
752

753
754
755
    /* Be clean */
    free(pid);
  }
756

757
  if (timer) TIMER_TOC(timer_do_ghost);
758
759
}

760
/**
761
 * @brief Unskip any tasks associated with active cells.
762
763
 *
 * @param c The cell.
764
 * @param e The engine.
765
 */
766
static void runner_do_unskip(struct cell *c, struct engine *e) {
767

768
769
770
  /* Ignore empty cells. */
  if (c->count == 0 && c->gcount == 0) return;

771
  /* Unskip any active tasks. */
772
  if (cell_is_active(c, e)) {
773
    const int forcerebuild = cell_unskip_tasks(c, &e->sched);
774
775
776
    if (forcerebuild) atomic_inc(&e->forcerebuild);
  }

777
  /* Recurse */
778
779
  if (c->split) {
    for (int k = 0; k < 8; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
780
      if (c->progeny[k] != NULL) {
Matthieu Schaller's avatar
Matthieu Schaller committed
781
        struct cell *cp = c->progeny[k];
782
        runner_do_unskip(cp, e);
783
784
785
      }
    }
  }
786
}
787

788
/**
789
 * @brief Mapper function to unskip active tasks.
790
791
792
793
794
 *
 * @param map_data An array of #cell%s.
 * @param num_elements Chunk size.
 * @param extra_data Pointer to an #engine.
 */
795
796
void runner_do_unskip_mapper(void *map_data, int num_elements,
                             void *extra_data) {
797

798
799
  struct engine *e = (struct engine *)extra_data;
  struct cell *cells = (struct cell *)map_data;
800

801
802
  for (int ind = 0; ind < num_elements; ind++) {
    struct cell *c = &cells[ind];
803
    if (c != NULL) runner_do_unskip(c, e);
804
  }
805
}
806
807
808
809
810
811
812
/**
 * @brief Drift particles in real space.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
813
void runner_do_drift_particles(struct runner *r, struct cell *c, int timer) {
814

815
  TIMER_TIC;
Matthieu Schaller's avatar
Matthieu Schaller committed
816

817
  cell_drift_particles(c, r->e);
818

819
  if (timer) TIMER_TOC(timer_drift);
820
}
821

822
/**
823
 * @brief Mapper function to drift ALL particle and g-particles forward in time.
824
825
826
827
828
829
830
 *
 * @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) {
831

832
833
834
835
836
  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];
837
    if (c != NULL && c->nodeID == e->nodeID) cell_drift_particles(c, e);
838
  }
839
}
840

841
842
843
844
845
846
847
/**
 * @brief Perform the first half-kick on all the active particles in a cell.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
848
void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
849

850
851
852
853
  const struct engine *e = r->e;
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
  struct gpart *restrict gparts = c->gparts;
854
  struct spart *restrict sparts = c->sparts;
855
856
  const int count = c->count;
  const int gcount = c->gcount;
857
  const int scount = c->scount;
858
  const integertime_t ti_current = e->ti_current;
859
  const double timeBase = e->timeBase;
860

861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
  TIMER_TIC;

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

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_kick1(r, c->progeny[k], 0);
  } else {

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

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

      /* If particle needs to be kicked */
      if (part_is_active(p, e)) {

        const integertime_t ti_step = get_integer_timestep(p->time_bin);
        const integertime_t ti_begin =
            get_integer_time_begin(ti_current, p->time_bin);

886
887
888
889
#ifdef SWIFT_DEBUG_CHECKS
        const integertime_t ti_end =
            get_integer_time_end(ti_current, p->time_bin);

890
        if (ti_end - ti_begin != ti_step)
891
892
893
894
          error(
              "Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
              "ti_step=%lld time_bin=%d ti_current=%lld",
              ti_end, ti_begin, ti_step, p->time_bin, ti_current);
895
896
#endif

897
        /* do the kick */
898
        kick_part(p, xp, ti_begin, ti_begin + ti_step / 2, timeBase);
899
900
901
      }
    }

902
    /* Loop over the gparts in this cell. */
903
904
905
906
907
908
    for (int k = 0; k < gcount; k++) {

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

      /* If the g-particle has no counterpart and needs to be kicked */
909
      if (gp->type == swift_type_dark_matter && gpart_is_active(gp, e)) {
910

911
912
913
914
        const integertime_t ti_step = get_integer_timestep(gp->time_bin);
        const integertime_t ti_begin =
            get_integer_time_begin(ti_current, gp->time_bin);

915
916
917
918
919
920
921
#ifdef SWIFT_DEBUG_CHECKS
        const integertime_t ti_end =
            get_integer_time_end(ti_current, gp->time_bin);

        if (ti_end - ti_begin != ti_step) error("Particle in wrong time-bin");
#endif

922
        /* do the kick */
923
        kick_gpart(gp, ti_begin, ti_begin + ti_step / 2, timeBase);
924
925
      }
    }
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950

    /* Loop over the star particles in this cell. */
    for (int k = 0; k < scount; k++) {

      /* Get a handle on the s-part. */
      struct spart *restrict sp = &sparts[k];

      /* If particle needs to be kicked */
      if (spart_is_active(sp, e)) {

        const integertime_t ti_step = get_integer_timestep(sp->time_bin);
        const integertime_t ti_begin =
            get_integer_time_begin(ti_current, sp->time_bin);

#ifdef SWIFT_DEBUG_CHECKS
        const integertime_t ti_end =
            get_integer_time_end(ti_current, sp->time_bin);

        if (ti_end - ti_begin != ti_step) error("Particle in wrong time-bin");
#endif

        /* do the kick */
        kick_spart(sp, ti_begin, ti_begin + ti_step / 2, timeBase);
      }
    }
951
  }
952

953
954
955
  if (timer) TIMER_TOC(timer_kick1);
}

956
957
958
/**
 * @brief Perform the second half-kick on all the active particles in a cell.
 *
959
 * Also prepares particles to be drifted.
960
961
962
963
964
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
965
void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
966
967
968

  const struct engine *e = r->e;
  const integertime_t ti_current = e->ti_current;
969
  const double timeBase = e->timeBase;
970
  const int count = c->count;
971
  const int gcount = c->gcount;
972
  const int scount = c->scount;
973
974
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
975
  struct gpart *restrict gparts = c->gparts;
976
  struct spart *restrict sparts = c->sparts;
977
978
979
980

  TIMER_TIC;

  /* Anything to do here? */
981
  if (!cell_is_active(c, e)) return;
982