runner.c 63.5 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 "chemistry.h"
46
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
47
#include "cooling.h"
48
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
49
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
50
#include "engine.h"
51
#include "error.h"
52
53
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
54
#include "hydro_properties.h"
55
#include "kick.h"
56
#include "minmax.h"
57
#include "runner_doiact_fft.h"
James Willis's avatar
James Willis committed
58
#include "runner_doiact_vec.h"
59
#include "scheduler.h"
60
#include "sort_part.h"
61
#include "sourceterms.h"
62
#include "space.h"
63
#include "stars.h"
64
65
#include "task.h"
#include "timers.h"
66
#include "timestep.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
67

68
69
70
71
#define TASK_LOOP_DENSITY 0
#define TASK_LOOP_GRADIENT 1
#define TASK_LOOP_FORCE 2
#define TASK_LOOP_LIMITER 3
72

73
/* Import the density loop functions. */
74
#define FUNCTION density
75
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
76
#include "runner_doiact.h"
77
78
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
79

80
/* Import the gradient loop functions (if required). */
81
82
#ifdef EXTRA_HYDRO_LOOP
#define FUNCTION gradient
83
#define FUNCTION_TASK_LOOP TASK_LOOP_GRADIENT
84
#include "runner_doiact.h"
85
86
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
87
88
#endif

89
/* Import the force loop functions. */
90
#define FUNCTION force
91
#define FUNCTION_TASK_LOOP TASK_LOOP_FORCE
92
#include "runner_doiact.h"
93
94
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
95

96
/* Import the gravity loop functions. */
97
#include "runner_doiact_fft.h"
98
#include "runner_doiact_grav.h"
99

Tom Theuns's avatar
Tom Theuns committed
100
/**
Tom Theuns's avatar
Tom Theuns committed
101
 * @brief Perform source terms
Tom Theuns's avatar
Tom Theuns committed
102
103
104
105
106
107
108
 *
 * @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;
109
  const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
Tom Theuns's avatar
Tom Theuns committed
110
  const double cell_width[3] = {c->width[0], c->width[1], c->width[2]};
Tom Theuns's avatar
Tom Theuns committed
111
  struct sourceterms *sourceterms = r->e->sourceterms;
112
  const int dimen = 3;
Tom Theuns's avatar
Tom Theuns committed
113
114
115
116
117
118
119

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

122
    if (count > 0) {
Tom Theuns's avatar
Tom Theuns committed
123

124
125
126
127
128
129
      /* 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
130
131
    }
  }
Tom Theuns's avatar
Tom Theuns committed
132
133
134
135

  if (timer) TIMER_TOC(timer_dosource);
}

Tom Theuns's avatar
Tom Theuns committed
136
137
138
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
139
140
141
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
142
 */
143
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
144

Matthieu Schaller's avatar
Matthieu Schaller committed
145
146
  struct gpart *restrict gparts = c->gparts;
  const int gcount = c->gcount;
147
148
149
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
150
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
151

152
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
153

154
  /* Anything to do here? */
155
  if (!cell_is_active_gravity(c, e)) return;
156

Tom Theuns's avatar
Tom Theuns committed
157
158
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
159
    for (int k = 0; k < 8; k++)
160
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
161
  } else {
162

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

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

169
      /* Is this part within the time step? */
170
      if (gpart_is_active(gp, e)) {
171
172
        external_gravity_acceleration(time, potential, constants, gp);
      }
173
    }
174
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
175

176
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
177
178
}

Stefan Arridge's avatar
Stefan Arridge committed
179
/**
180
181
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
182
183
184
185
186
187
188
 *
 * @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) {

189
  const struct engine *e = r->e;
190
191
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
192
193
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
194
  const struct unit_system *us = e->internal_units;
195
  const double time_base = e->time_base;
196
197
198
199
  const integertime_t ti_current = e->ti_current;
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
  const int count = c->count;
Stefan Arridge's avatar
Stefan Arridge committed
200
201
202

  TIMER_TIC;

203
  /* Anything to do here? */
204
  if (!cell_is_active_hydro(c, e)) return;
205

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

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

215
216
217
      /* 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
218

219
      if (part_is_active(p, e)) {
220

221
222
223
224
225
226
227
228
229
230
        double dt_cool;
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(ti_current + 1, p->time_bin);
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
        }
231

232
        /* Let's cool ! */
233
        cooling_cool_part(constants, us, cosmo, cooling_func, p, xp, dt_cool);
234
      }
Stefan Arridge's avatar
Stefan Arridge committed
235
236
237
238
239
240
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
319
320
321
322
323
324
325
326
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
 * Debugging function.
 *
 * @param c The #cell to check.
 * @param flags The sorting flags to check.
 */
327
void runner_check_sorts(struct cell *c, int flags) {
Matthieu Schaller's avatar
Matthieu Schaller committed
328
329

#ifdef SWIFT_DEBUG_CHECKS
Pedro Gonnet's avatar
Pedro Gonnet committed
330
  if (flags & ~c->sorted) error("Inconsistent sort flags (downward)!");
331
332
  if (c->split)
    for (int k = 0; k < 8; k++)
333
      if (c->progeny[k] != NULL && c->progeny[k]->count > 0)
334
        runner_check_sorts(c->progeny[k], c->sorted);
Matthieu Schaller's avatar
Matthieu Schaller committed
335
336
337
#else
  error("Calling debugging code without debugging flag activated.");
#endif
338
339
}

Pedro Gonnet's avatar
Pedro Gonnet committed
340
341
342
343
344
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
345
 * @param flags Cell flag.
346
347
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
348
349
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
350
 */
351
352
void runner_do_sort(struct runner *r, struct cell *c, int flags, int cleanup,
                    int clock) {
353
354

  struct entry *fingers[8];
355
  const int count = c->count;
356
357
  const struct part *parts = c->parts;
  struct xpart *xparts = c->xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
358
  float buff[8];
359

360
361
  TIMER_TIC;

362
363
  /* We need to do the local sorts plus whatever was requested further up. */
  flags |= c->do_sort;
364
365
366
367
368
  if (cleanup) {
    c->sorted = 0;
  } else {
    flags &= ~c->sorted;
  }
369
  if (flags == 0 && !c->do_sub_sort) return;
370
371

  /* Check that the particles have been moved to the current time */
Pedro Gonnet's avatar
Pedro Gonnet committed
372
  if (flags && !cell_are_part_drifted(c, r->e))
373
    error("Sorting un-drifted cell c->nodeID=%d", c->nodeID);
Pedro Gonnet's avatar
Pedro Gonnet committed
374

375
376
377
378
379
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
  runner_check_sorts(c, c->sorted);

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
380
381
382
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
    if (finger->sorted & ~c->sorted) error("Inconsistent sort flags (upward).");
383
  }
384

385
386
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
387
  if (c->sorted == 0) c->ti_sort = r->e->ti_current;
388
#endif
389

390
391
392
393
394
395
396
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
    if ((flags & (1 << j)) && c->sort[j] == NULL) {
      if ((c->sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                               (count + 1))) == NULL)
        error("Failed to allocate sort memory.");
    }
397
398
399
400
401
402
  }

  /* Does this cell have any progeny? */
  if (c->split) {

    /* Fill in the gaps within the progeny. */
403
    float dx_max_sort = 0.0f;
404
    float dx_max_sort_old = 0.0f;
405
    for (int k = 0; k < 8; k++) {
406
      if (c->progeny[k] != NULL && c->progeny[k]->count > 0) {
407
408
409
410
411
        /* Only propagate cleanup if the progeny is stale. */
        runner_do_sort(r, c->progeny[k], flags,
                       cleanup && (c->progeny[k]->dx_max_sort >
                                   space_maxreldx * c->progeny[k]->dmin),
                       0);
412
        dx_max_sort = max(dx_max_sort, c->progeny[k]->dx_max_sort);
413
        dx_max_sort_old = max(dx_max_sort_old, c->progeny[k]->dx_max_sort_old);
414
      }
415
    }
416
    c->dx_max_sort = dx_max_sort;
417
    c->dx_max_sort_old = dx_max_sort_old;
418
419

    /* Loop over the 13 different sort arrays. */
420
    for (int j = 0; j < 13; j++) {
421
422
423
424
425

      /* Has this sort array been flagged? */
      if (!(flags & (1 << j))) continue;

      /* Init the particle index offsets. */
426
      int off[8];
427
428
      off[0] = 0;
      for (int k = 1; k < 8; k++)
429
430
431
432
433
434
        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. */
435
      int inds[8];
436
      for (int k = 0; k < 8; k++) {
437
438
        inds[k] = k;
        if (c->progeny[k] != NULL && c->progeny[k]->count > 0) {
439
          fingers[k] = c->progeny[k]->sort[j];
440
441
442
443
444
445
446
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
447
448
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
449
          if (buff[inds[k]] < buff[inds[i]]) {
450
            int temp_i = inds[i];
451
452
453
454
455
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
456
      struct entry *finger = c->sort[j];
457
      for (int ind = 0; ind < count; ind++) {
458
459
460
461
462
463
464
465
466
467

        /* 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. */
468
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
469
          int temp_i = inds[k - 1];
470
471
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
472
        }
473

474
475
476
      } /* Merge. */

      /* Add a sentinel. */
477
478
      c->sort[j][count].d = FLT_MAX;
      c->sort[j][count].i = 0;
479
480

      /* Mark as sorted. */
481
      atomic_or(&c->sorted, 1 << j);
482
483
484
485
486
487
488
489

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

490
    /* Reset the sort distance */
491
    if (c->sorted == 0) {
492
493
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
494
        error("Have non-NULL xparts in foreign cell");
495
#endif
496
497
498
499
500
501
502
503

      /* And the individual sort distances if we are a local cell */
      if (xparts != NULL) {
        for (int k = 0; k < count; k++) {
          xparts[k].x_diff_sort[0] = 0.0f;
          xparts[k].x_diff_sort[1] = 0.0f;
          xparts[k].x_diff_sort[2] = 0.0f;
        }
504
      }
505
506
      c->dx_max_sort_old = 0.f;
      c->dx_max_sort = 0.f;
507
508
    }

509
    /* Fill the sort array. */
510
    for (int k = 0; k < count; k++) {
511
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
512
      for (int j = 0; j < 13; j++)
513
        if (flags & (1 << j)) {
514
515
516
517
          c->sort[j][k].i = k;
          c->sort[j][k].d = px[0] * runner_shift[j][0] +
                            px[1] * runner_shift[j][1] +
                            px[2] * runner_shift[j][2];
518
        }
519
    }
520
521

    /* Add the sentinel and sort. */
522
    for (int j = 0; j < 13; j++)
523
      if (flags & (1 << j)) {
524
525
526
        c->sort[j][count].d = FLT_MAX;
        c->sort[j][count].i = 0;
        runner_do_sort_ascending(c->sort[j], count);
527
        atomic_or(&c->sorted, 1 << j);
528
529
530
      }
  }

531
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
532
  /* Verify the sorting. */
533
  for (int j = 0; j < 13; j++) {
534
    if (!(flags & (1 << j))) continue;
535
    struct entry *finger = c->sort[j];
536
    for (int k = 1; k < count; k++) {
537
538
539
540
541
      if (finger[k].d < finger[k - 1].d)
        error("Sorting failed, ascending array.");
      if (finger[k].i >= count) error("Sorting failed, indices borked.");
    }
  }
Pedro Gonnet's avatar
Pedro Gonnet committed
542

543
544
545
546
  /* Make sure the sort flags are consistent (downward). */
  runner_check_sorts(c, flags);

  /* Make sure the sort flags are consistent (upward). */
Pedro Gonnet's avatar
Pedro Gonnet committed
547
548
549
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
    if (finger->sorted & ~c->sorted) error("Inconsistent sort flags.");
550
  }
551
#endif
552

553
554
555
556
557
  /* Clear the cell's sort flags. */
  c->do_sort = 0;
  c->do_sub_sort = 0;
  c->requires_sorts = 0;

558
559
560
  if (clock) TIMER_TOC(timer_dosort);
}

561
/**
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
 * @brief Initialize the multipoles before the gravity calculation.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer 1 if the time is to be recorded.
 */
void runner_do_init_grav(struct runner *r, struct cell *c, int timer) {

  const struct engine *e = r->e;

  TIMER_TIC;

#ifdef SWIFT_DEBUG_CHECKS
  if (!(e->policy & engine_policy_self_gravity))
    error("Grav-init task called outside of self-gravity calculation");
#endif

  /* Anything to do here? */
580
  if (!cell_is_active_gravity(c, e)) return;
581
582

  /* Reset the gravity acceleration tensors */
583
  gravity_field_tensors_init(&c->multipole->pot, e->ti_current);
584
585
586
587
588
589
590
591
592
593
594

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

  if (timer) TIMER_TOC(timer_init_grav);
}

595
/**
596
597
598
599
600
 * @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.
601
 * @param timer Are we timing this ?
602
 */
603
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
604

605
#ifdef EXTRA_HYDRO_LOOP
606

607
608
  struct part *restrict parts = c->parts;
  const int count = c->count;
609
  const struct engine *e = r->e;
610

611
612
  TIMER_TIC;

613
  /* Anything to do here? */
614
  if (!cell_is_active_hydro(c, e)) return;
615

616
617
618
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
619
      if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k], 0);
620
621
622
623
624
625
626
627
  } 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];

628
      if (part_is_active(p, e)) {
629
630
631
632
633
634

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

636
637
  if (timer) TIMER_TOC(timer_do_extra_ghost);

638
639
#else
  error("SWIFT was not compiled with the extra hydro loop activated.");
640
#endif
641
}
642

643
/**
644
645
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
646
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
647
 * @param r The runner thread.
648
 * @param c The cell.
649
 * @param timer Are we timing this ?
650
 */
651
void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
652

Matthieu Schaller's avatar
Matthieu Schaller committed
653
654
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
655
  const struct engine *e = r->e;
656
  const struct space *s = e->s;
657
658
  const struct hydro_space *hs = &s->hs;
  const struct cosmology *cosmo = e->cosmology;
659
  const struct chemistry_data *chemistry = e->chemistry;
660
  const float hydro_h_max = e->hydro_properties->h_max;
661
662
663
  const float eps = e->hydro_properties->h_tolerance;
  const float hydro_eta_dim =
      pow_dimension(e->hydro_properties->eta_neighbours);
664
  const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
665
  int redo = 0, count = 0;
666

667
668
  TIMER_TIC;

669
  /* Anything to do here? */
670
  if (!cell_is_active_hydro(c, e)) return;
671

672
673
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
674
    for (int k = 0; k < 8; k++)
675
676
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k], 0);
  } else {
677

678
    /* Init the list of active particles that have to be updated. */
679
    int *pid = NULL;
680
    if ((pid = (int *)malloc(sizeof(int) * c->count)) == NULL)
681
      error("Can't allocate memory for pid.");
682
683
684
685
686
    for (int k = 0; k < c->count; k++)
      if (part_is_active(&parts[k], e)) {
        pid[count] = k;
        ++count;
      }
687

688
689
690
    /* While there are particles that need to be updated... */
    for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
         num_reruns++) {
691

692
693
      /* Reset the redo-count. */
      redo = 0;
694

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

698
        /* Get a direct pointer on the part. */
699
700
        struct part *p = &parts[pid[i]];
        struct xpart *xp = &xparts[pid[i]];
701

702
#ifdef SWIFT_DEBUG_CHECKS
703
        /* Is this part within the timestep? */
704
705
706
        if (!part_is_active(p, e)) error("Ghost applied to inactive particle");
#endif

707
708
709
710
711
        /* Get some useful values */
        const float h_old = p->h;
        const float h_old_dim = pow_dimension(h_old);
        const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
        float h_new;
712

713
        if (p->density.wcount == 0.f) { /* No neighbours case */
714

715
716
717
          /* Double h and try again */
          h_new = 2.f * h_old;
        } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
718

719
          /* Finish the density calculation */
720
          hydro_end_density(p, cosmo);
721
          chemistry_end_density(p, chemistry, cosmo);
722

723
724
725
726
727
728
729
          /* Compute one step of the Newton-Raphson scheme */
          const float n_sum = p->density.wcount * h_old_dim;
          const float n_target = hydro_eta_dim;
          const float f = n_sum - n_target;
          const float f_prime =
              p->density.wcount_dh * h_old_dim +
              hydro_dimension * p->density.wcount * h_old_dim_minus_one;
730

731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
          h_new = h_old - f / f_prime;

#ifdef SWIFT_DEBUG_CHECKS
          if ((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old))
            error(
                "Smoothing length correction not going in the right direction");
#endif

          /* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
          h_new = min(h_new, 2.f * h_old);
          h_new = max(h_new, 0.5f * h_old);
        }

        /* Check whether the particle has an inappropriate smoothing length */
        if (fabsf(h_new - h_old) > eps * h_old) {
746

747
          /* Ok, correct then */
748
          p->h = h_new;
749

750
751
          /* If below the absolute maximum, try again */
          if (p->h < hydro_h_max) {
752

753
754
755
            /* Flag for another round of fun */
            pid[redo] = pid[i];
            redo += 1;
756

757
            /* Re-initialise everything */
758
            hydro_init_part(p, hs);
759
            chemistry_init_part(p, chemistry);
760
761
762

            /* Off we go ! */
            continue;
763

764
765
766
767
          } else {

            /* Ok, this particle is a lost cause... */
            p->h = hydro_h_max;
768
769
770

            /* Do some damage control if no neighbours at all were found */
            if (p->density.wcount == kernel_root * kernel_norm)
771
              hydro_part_has_no_neighbours(p, xp, cosmo);
772
          }
773
        }
774

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

777
        /* As of here, particle force variables will be set. */
778

779
        /* Compute variables required for the force loop */
780
        hydro_prepare_force(p, xp, cosmo);
781

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

785
786
        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
787
788
      }

789
790
      /* We now need to treat the particles whose smoothing length had not
       * converged again */
791

792
793
794
      /* Re-set the counter for the next loop (potentially). */
      count = redo;
      if (count > 0) {
795

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

799
800
          /* 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
801

802
803
804
805
#ifdef SWIFT_DEBUG_CHECKS
            if (l->t->ti_run < r->e->ti_current)
              error("Density task should have been run.");
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
806

807
808
            /* Self-interaction? */
            if (l->t->type == task_type_self)
809
              runner_doself_subset_branch_density(r, finger, parts, pid, count);
Matthieu Schaller's avatar
Matthieu Schaller committed
810

811
812
            /* Otherwise, pair interaction? */
            else if (l->t->type == task_type_pair) {
813

814
815
              /* Left or right? */
              if (l->t->ci == finger)
James Willis's avatar
James Willis committed
816
817
                runner_dopair_subset_branch_density(r, finger, parts, pid,
                                                    count, l->t->cj);
818
              else
James Willis's avatar
James Willis committed
819
820
                runner_dopair_subset_branch_density(r, finger, parts, pid,
                                                    count, l->t->ci);
821
            }
822

823
824
825
826
            /* 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);
827

828
829
830
831
832
833
834
835
836
837
838
            /* 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);
            }
839
840
841
          }
        }
      }
842
    }
843

844
845
#ifdef SWIFT_DEBUG_CHECKS
    if (count) {
846
      error("Smoothing length failed to converge on %i particles.", count);
847
848
    }
#else
849
    if (count)
850
      error("Smoothing length failed to converge on %i particles.", count);
851
#endif
852

853
854
855
    /* Be clean */
    free(pid);
  }
856

857
  if (timer) TIMER_TOC(timer_do_ghost);
858
859
}

860
/**
861
 * @brief Unskip any hydro tasks associated with active cells.
862
863
 *
 * @param c The cell.
864
 * @param e The engine.
865
 */
866
static void runner_do_unskip_hydro(struct cell *c, struct engine *e) {
867

868
  /* Ignore empty cells. */
869
  if (c->count == 0) return;
870

871
  /* Skip inactive cells. */
872
  if (!cell_is_active_hydro(c, e)) return;
873

874
  /* Recurse */
875
876
  if (c->split) {
    for (int k = 0; k < 8; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
877
      if (c->progeny[k] != NULL) {
Matthieu Schaller's avatar
Matthieu Schaller committed
878
        struct cell *cp = c->progeny[k];
879
        runner_do_unskip_hydro(cp, e);
880
881
882
      }
    }
  }