runner.c 40.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"
James Willis's avatar
James Willis committed
56
#include "runner_doiact_vec.h"
57
#include "scheduler.h"
58
#include "sourceterms.h"
59
60
61
#include "space.h"
#include "task.h"
#include "timers.h"
62
#include "timestep.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
  const int count = c->count;
209
  // 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
      /* Kick has already updated ti_end, so need to check ti_begin */
231
232
      // if (p->ti_begin == ti_current) {
      {
233

234
235
        // const double dt = (p->ti_end - p->ti_begin) * timeBase;
        const double dt = 1. * timeBase;  // MATTHIEU
236

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

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

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

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

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

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

    /* 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
419
        }
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
      } /* 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
445
446
447
          sort[j * (count + 1) + k].d = px[0] * runner_shift[j][0] +
                                        px[1] * runner_shift[j][1] +
                                        px[2] * runner_shift[j][2];
448
        }
449
    }
450
451
452
453
454
455

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

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

  if (clock) TIMER_TOC(timer_dosort);
}

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

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

  TIMER_TIC;
493

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

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

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

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

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

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

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

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

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

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

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

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

543
#ifdef EXTRA_HYDRO_LOOP
544

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

549
550
  TIMER_TIC;

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

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

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

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

574
575
  if (timer) TIMER_TOC(timer_do_extra_ghost);

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
591
592
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
593
  int redo, count = c->count;
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

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
          /* Compute variables required for the force loop */
Matthieu Schaller's avatar
Matthieu Schaller committed
676
          hydro_prepare_force(p, xp);
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
/**
748
 * @brief Unskip any tasks associated with active cells.
749
750
 *
 * @param c The cell.
751
 * @param e The engine.
752
 */
753
static void runner_do_unskip(struct cell *c, struct engine *e) {
754

755
  /* Unskip any active tasks. */
756
  if (cell_is_active(c, e)) {
757
    const int forcerebuild = cell_unskip_tasks(c, &e->sched);
758
759
760
    if (forcerebuild) atomic_inc(&e->forcerebuild);
  }

761
  /* Recurse */
762
763
  if (c->split) {
    for (int k = 0; k < 8; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
764
      if (c->progeny[k] != NULL) {
Matthieu Schaller's avatar
Matthieu Schaller committed
765
        struct cell *cp = c->progeny[k];
766
        runner_do_unskip(cp, e);
767
768
769
      }
    }
  }
770
}
771

772
/**
773
 * @brief Mapper function to unskip active tasks.
774
775
776
777
778
 *
 * @param map_data An array of #cell%s.
 * @param num_elements Chunk size.
 * @param extra_data Pointer to an #engine.
 */
779
780
void runner_do_unskip_mapper(void *map_data, int num_elements,
                             void *extra_data) {
781

782
783
  struct engine *e = (struct engine *)extra_data;
  struct cell *cells = (struct cell *)map_data;
784

785
786
  for (int ind = 0; ind < num_elements; ind++) {
    struct cell *c = &cells[ind];
787
    if (c != NULL) runner_do_unskip(c, e);
788
  }
789
}
790
791
792
793
794
795
796
/**
 * @brief Drift particles in real space.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
797
void runner_do_drift(struct runner *r, struct cell *c, int timer) {
798

799
  TIMER_TIC;
Matthieu Schaller's avatar
Matthieu Schaller committed
800

801
  cell_drift(c, r->e);
802

803
  if (timer) TIMER_TOC(timer_drift);
804
}
805

806
/**
807
 * @brief Mapper function to drift ALL particle and g-particles forward in time.
808
809
810
811
812
813
814
 *
 * @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) {
815

816
817
818
819
820
  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];
821
    if (c != NULL && c->nodeID == e->nodeID) cell_drift(c, e);
822
  }
823
}
824

825
826
827
828
829
830
831
/**
 * @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 ?
 */
832
void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
833

834
835
836
837
838
839
840
  const struct engine *e = r->e;
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
  struct gpart *restrict gparts = c->gparts;
  const int count = c->count;
  const int gcount = c->gcount;
  const integertime_t ti_current = e->ti_current;
841
  const double timeBase = e->timeBase;
842

843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
  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);

868
869
870
871
872
873
874
#ifdef SWIFT_DEBUG_CHECKS
        const integertime_t ti_end =
            get_integer_time_end(ti_current, p->time_bin);

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

875
        /* do the kick */
876
        kick_part(p, xp, ti_begin, ti_begin + ti_step / 2, timeBase);
877
878
879
      }
    }

880
    /* Loop over the gparts in this cell. */
881
882
883
884
885
886
887
888
    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 */
      if (gp->id_or_neg_offset > 0 && gpart_is_active(gp, e)) {

889
890
891
892
        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);

893
894
895
896
897
898
899
#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

900
        /* do the kick */
901
        kick_gpart(gp, ti_begin, ti_begin + ti_step / 2, timeBase);
902
903
904
905
906
907
      }
    }
  }
  if (timer) TIMER_TOC(timer_kick1);
}

908
909
910
/**
 * @brief Perform the second half-kick on all the active particles in a cell.
 *
911
 * Also prepares particles to be drifted.
912
913
914
915
916
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
917
void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
918
919
920

  const struct engine *e = r->e;
  const integertime_t ti_current = e->ti_current;
921
  const double timeBase = e->timeBase;
922
  const int count = c->count;
923
  const int gcount = c->gcount;
924
925
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
926
  struct gpart *restrict gparts = c->gparts;
927
928
929
930
931
  const double const_G = e->physical_constants->const_newton_G;

  TIMER_TIC;

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

934
935
936
937
938
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_kick2(r, c->progeny[k], 0);
  } else {
939

940
    /* Loop over the particles in this cell. */
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
    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)) {

        /* First, finish the force loop */
        hydro_end_force(p);
        if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);

        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);

958
959
960
961
962
#ifdef SWIFT_DEBUG_CHECKS
        if (ti_begin + ti_step != ti_current)
          error("Particle in wrong time-bin");
#endif

963
        /* Finish the time-step with a second half-kick */
964
        kick_part(p, xp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
965
966
967

        /* Prepare the values to be drifted */
        hydro_reset_predicted_values(p, xp);
968
969
970
971
972
973
974
975
976
      }
    }

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

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

977
978
      /* If the g-particle has no counterpart and needs to be kicked */
      if (gp->id_or_neg_offset > 0 && gpart_is_active(gp, e)) {
979

980
981
        /* First, finish the force loop */
        gravity_end_force(gp, const_G);
982

983
984
985
        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);
986

987
#ifdef SWIFT_DEBUG_CHECKS
988
989
        if (ti_begin + ti_step != ti_current)
          error("Particle in wrong time-bin");