runner.c 75.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"
James Willis's avatar
James Willis committed
57
#include "runner_doiact_vec.h"
58
#include "scheduler.h"
59
#include "sort_part.h"
60
#include "sourceterms.h"
61
#include "space.h"
62
#include "space_getsid.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_grav.h"
98

Loic Hausammann's avatar
Loic Hausammann committed
99
100
/* Import the stars loop functions. */
#include "runner_doiact_stars.h"
101

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

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

124
    if (count > 0) {
Tom Theuns's avatar
Tom Theuns committed
125

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

  if (timer) TIMER_TOC(timer_dosource);
}

138
139
140
141
142
143
144
145
/**
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
Loic Hausammann's avatar
Loic Hausammann committed
146
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
147
148
149
150

  struct spart *restrict sparts = c->sparts;
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
Loic Hausammann's avatar
Loic Hausammann committed
151
  const struct stars_props *stars_properties = e->stars_properties;
Loic Hausammann's avatar
Loic Hausammann committed
152
  const float stars_h_max = stars_properties->h_max;
153
  const float eps = stars_properties->h_tolerance;
Loic Hausammann's avatar
Loic Hausammann committed
154
  const float stars_eta_dim =
155
156
157
      pow_dimension(stars_properties->eta_neighbours);
  const int max_smoothing_iter = stars_properties->max_smoothing_iterations;
  int redo = 0, scount = 0;
158
159
160
161

  TIMER_TIC;

  /* Anything to do here? */
Loic Hausammann's avatar
Loic Hausammann committed
162
  if (!cell_is_active_stars(c, e)) return;
163
164
165
166

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
Loic Hausammann's avatar
Loic Hausammann committed
167
      if (c->progeny[k] != NULL) runner_do_stars_ghost(r, c->progeny[k], 0);
168
169
170
171
172
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
    if ((sid = (int *)malloc(sizeof(int) * c->scount)) == NULL)
Loic Hausammann's avatar
Loic Hausammann committed
173
      error("Can't allocate memory for sid.");
174
175
    for (int k = 0; k < c->scount; k++)
      if (spart_is_active(&sparts[k], e)) {
176
177
        sid[scount] = k;
        ++scount;
178
179
180
      }

    /* While there are particles that need to be updated... */
181
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
182
183
184
185
186
187
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
188
      for (int i = 0; i < scount; i++) {
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204

        /* Get a direct pointer on the part. */
        struct spart *sp = &sparts[sid[i]];

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
        if (!spart_is_active(sp, e)) error("Ghost applied to inactive particle");
#endif

        /* Get some useful values */
        const float h_old = sp->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;
        int has_no_neighbours = 0;

205
        if (sp->density.wcount == 0.f) { /* No neighbours case */
206
207
208
209
210
211
212
213
214

          /* Flag that there were no neighbours */
          has_no_neighbours = 1;

          /* Double h and try again */
          h_new = 2.f * h_old;
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
215
          stars_end_density(sp, cosmo);
216
217

          /* Compute one step of the Newton-Raphson scheme */
218
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
219
          const float n_target = stars_eta_dim;
220
221
          const float f = n_sum - n_target;
          const float f_prime =
222
223
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244

          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
#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) {

          /* Ok, correct then */
          sp->h = h_new;

          /* If below the absolute maximum, try again */
Loic Hausammann's avatar
Loic Hausammann committed
245
          if (sp->h < stars_h_max) {
246
247
248
249
250
251

            /* Flag for another round of fun */
            sid[redo] = sid[i];
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
252
            stars_init_spart(sp);
253
254
255
256
257
258
259

            /* Off we go ! */
            continue;

          } else {

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
260
            sp->h = stars_h_max;
261
262
263

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
264
              stars_spart_has_no_neighbours(sp, cosmo);
265
266
267
268
            }
          }
        }

269
	/* We now have a particle whose smoothing length has converged */
270

Loic Hausammann's avatar
Loic Hausammann committed
271
        /* Compute the stellar evolution  */
272
        stars_evolve_spart(sp, stars_properties, cosmo);
273
274
275
276
277
278
279

      }

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

      /* Re-set the counter for the next loop (potentially). */
280
281
      scount = redo;
      if (scount > 0) {
282
283
284
285
286
287
288
289
290
291
292
293
294
295

        /* Climb up the cell hierarchy. */
        for (struct cell *finger = c; finger != NULL; finger = finger->parent) {

          /* Run through this cell's density interactions. */
          for (struct link *l = finger->density; l != NULL; l = l->next) {

#ifdef SWIFT_DEBUG_CHECKS
            if (l->t->ti_run < r->e->ti_current)
              error("Density task should have been run.");
#endif

            /* Self-interaction? */
            if (l->t->type == task_type_self)
Loic Hausammann's avatar
Loic Hausammann committed
296
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid, scount);
297
298
299
300
301
302

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

              /* Left or right? */
              if (l->t->ci == finger)
Loic Hausammann's avatar
Loic Hausammann committed
303
                runner_dopair_subset_branch_stars_density(r, finger, sparts, sid,
304
                                                    scount, l->t->cj);
305
              else
Loic Hausammann's avatar
Loic Hausammann committed
306
                runner_dopair_subset_branch_stars_density(r, finger, sparts, sid,
307
                                                    scount, l->t->ci);
308
309
310
311
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
Loic Hausammann's avatar
Loic Hausammann committed
312
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount, NULL,
313
314
315
316
317
318
319
                                          -1, 1);

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

              /* Left or right? */
              if (l->t->ci == finger)
Loic Hausammann's avatar
Loic Hausammann committed
320
                runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
321
322
                                            l->t->cj, -1, 1);
              else
Loic Hausammann's avatar
Loic Hausammann committed
323
                runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
324
325
326
327
328
329
330
                                            l->t->ci, -1, 1);
            }
          }
        }
      }
    }

331
332
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
333
334
335
336
337
338
    }

    /* Be clean */
    free(sid);
  }

Loic Hausammann's avatar
Loic Hausammann committed
339
  if (timer) TIMER_TOC(timer_do_stars_ghost);
340
341
342
}


Tom Theuns's avatar
Tom Theuns committed
343
344
345
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
346
347
348
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
349
 */
350
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
351

Matthieu Schaller's avatar
Matthieu Schaller committed
352
353
  struct gpart *restrict gparts = c->gparts;
  const int gcount = c->gcount;
354
355
356
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
357
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
358

359
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
360

361
  /* Anything to do here? */
362
  if (!cell_is_active_gravity(c, e)) return;
363

Tom Theuns's avatar
Tom Theuns committed
364
365
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
366
    for (int k = 0; k < 8; k++)
367
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
368
  } else {
369

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

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

376
      /* Is this part within the time step? */
377
      if (gpart_is_active(gp, e)) {
378
379
        external_gravity_acceleration(time, potential, constants, gp);
      }
380
    }
381
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
382

383
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
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
420
/**
 * @brief Calculate gravity accelerations from the periodic mesh
 *
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
 */
void runner_do_grav_mesh(struct runner *r, struct cell *c, int timer) {

  struct gpart *restrict gparts = c->gparts;
  const int gcount = c->gcount;
  const struct engine *e = r->e;

#ifdef SWIFT_DEBUG_CHECKS
  if (!e->s->periodic) error("Calling mesh forces in non-periodic mode.");
#endif

  TIMER_TIC;

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

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

    /* Get the forces from the gravity mesh */
    pm_mesh_interpolate_forces(e->mesh, e, gparts, gcount);
  }

  if (timer) TIMER_TOC(timer_dograv_mesh);
}

Stefan Arridge's avatar
Stefan Arridge committed
421
/**
422
423
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
424
425
426
427
428
429
430
 *
 * @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) {

431
  const struct engine *e = r->e;
432
433
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
434
435
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
436
  const struct unit_system *us = e->internal_units;
437
  const double time_base = e->time_base;
438
439
440
441
  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
442
443
444

  TIMER_TIC;

445
  /* Anything to do here? */
446
  if (!cell_is_active_hydro(c, e)) return;
447

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

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

457
458
459
      /* 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
460

461
      if (part_is_active(p, e)) {
462

463
464
465
466
467
468
469
470
471
472
        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);
        }
473

474
        /* Let's cool ! */
475
        cooling_cool_part(constants, us, cosmo, cooling_func, p, xp, dt_cool);
476
      }
Stefan Arridge's avatar
Stefan Arridge committed
477
478
479
480
481
482
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
483
484
485
486
487
488
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
489
void runner_do_sort_ascending(struct entry *sort, int N) {
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543

  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
544
        }
545
546
547
548
549
550
551
552
553
554
555
556
      } 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
557
    }
558
559
560
  }
}

Matthieu Schaller's avatar
Matthieu Schaller committed
561
562
563
564
565
566
567
568
/**
 * @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.
 */
569
void runner_check_sorts(struct cell *c, int flags) {
Matthieu Schaller's avatar
Matthieu Schaller committed
570
571

#ifdef SWIFT_DEBUG_CHECKS
Pedro Gonnet's avatar
Pedro Gonnet committed
572
  if (flags & ~c->sorted) error("Inconsistent sort flags (downward)!");
573
574
  if (c->split)
    for (int k = 0; k < 8; k++)
575
      if (c->progeny[k] != NULL && c->progeny[k]->count > 0)
576
        runner_check_sorts(c->progeny[k], c->sorted);
Matthieu Schaller's avatar
Matthieu Schaller committed
577
578
579
#else
  error("Calling debugging code without debugging flag activated.");
#endif
580
581
}

Pedro Gonnet's avatar
Pedro Gonnet committed
582
583
584
585
586
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
587
 * @param flags Cell flag.
588
589
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
590
591
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
592
 */
593
594
void runner_do_sort(struct runner *r, struct cell *c, int flags, int cleanup,
                    int clock) {
595
596

  struct entry *fingers[8];
597
  const int count = c->count;
598
599
  const struct part *parts = c->parts;
  struct xpart *xparts = c->xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
600
  float buff[8];
601

602
603
  TIMER_TIC;

604
605
  /* We need to do the local sorts plus whatever was requested further up. */
  flags |= c->do_sort;
606
607
608
609
610
  if (cleanup) {
    c->sorted = 0;
  } else {
    flags &= ~c->sorted;
  }
611
  if (flags == 0 && !c->do_sub_sort) return;
612
613

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

617
618
619
620
621
#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
622
623
624
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
    if (finger->sorted & ~c->sorted) error("Inconsistent sort flags (upward).");
625
  }
626

627
628
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
629
  if (c->sorted == 0) c->ti_sort = r->e->ti_current;
630
#endif
631

632
633
634
635
636
637
638
  /* 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.");
    }
639
640
641
642
643
644
  }

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

    /* Fill in the gaps within the progeny. */
645
    float dx_max_sort = 0.0f;
646
    float dx_max_sort_old = 0.0f;
647
    for (int k = 0; k < 8; k++) {
648
      if (c->progeny[k] != NULL && c->progeny[k]->count > 0) {
649
650
651
652
653
        /* 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);
654
        dx_max_sort = max(dx_max_sort, c->progeny[k]->dx_max_sort);
655
        dx_max_sort_old = max(dx_max_sort_old, c->progeny[k]->dx_max_sort_old);
656
      }
657
    }
658
    c->dx_max_sort = dx_max_sort;
659
    c->dx_max_sort_old = dx_max_sort_old;
660
661

    /* Loop over the 13 different sort arrays. */
662
    for (int j = 0; j < 13; j++) {
663
664
665
666
667

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

      /* Init the particle index offsets. */
668
      int off[8];
669
670
      off[0] = 0;
      for (int k = 1; k < 8; k++)
671
672
673
674
675
676
        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. */
677
      int inds[8];
678
      for (int k = 0; k < 8; k++) {
679
680
        inds[k] = k;
        if (c->progeny[k] != NULL && c->progeny[k]->count > 0) {
681
          fingers[k] = c->progeny[k]->sort[j];
682
683
684
685
686
687
688
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
689
690
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
691
          if (buff[inds[k]] < buff[inds[i]]) {
692
            int temp_i = inds[i];
693
694
695
696
697
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
698
      struct entry *finger = c->sort[j];
699
      for (int ind = 0; ind < count; ind++) {
700
701
702
703
704
705
706
707
708
709

        /* 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. */
710
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
711
          int temp_i = inds[k - 1];
712
713
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
714
        }
715

716
717
718
      } /* Merge. */

      /* Add a sentinel. */
719
720
      c->sort[j][count].d = FLT_MAX;
      c->sort[j][count].i = 0;
721
722

      /* Mark as sorted. */
723
      atomic_or(&c->sorted, 1 << j);
724
725
726
727
728
729
730
731

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

732
    /* Reset the sort distance */
733
    if (c->sorted == 0) {
734
735
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
736
        error("Have non-NULL xparts in foreign cell");
737
#endif
738
739
740
741
742
743
744
745

      /* 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;
        }
746
      }
747
748
      c->dx_max_sort_old = 0.f;
      c->dx_max_sort = 0.f;
749
750
    }

751
    /* Fill the sort array. */
752
    for (int k = 0; k < count; k++) {
753
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
754
      for (int j = 0; j < 13; j++)
755
        if (flags & (1 << j)) {
756
757
758
759
          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];
760
        }
761
    }
762
763

    /* Add the sentinel and sort. */
764
    for (int j = 0; j < 13; j++)
765
      if (flags & (1 << j)) {
766
767
768
        c->sort[j][count].d = FLT_MAX;
        c->sort[j][count].i = 0;
        runner_do_sort_ascending(c->sort[j], count);
769
        atomic_or(&c->sorted, 1 << j);
770
771
772
      }
  }

773
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
774
  /* Verify the sorting. */
775
  for (int j = 0; j < 13; j++) {
776
    if (!(flags & (1 << j))) continue;
777
    struct entry *finger = c->sort[j];
778
    for (int k = 1; k < count; k++) {
779
780
781
782
783
      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
784

785
786
787
788
  /* 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
789
790
791
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
    if (finger->sorted & ~c->sorted) error("Inconsistent sort flags.");
792
  }
793
#endif
794

795
796
797
798
799
  /* Clear the cell's sort flags. */
  c->do_sort = 0;
  c->do_sub_sort = 0;
  c->requires_sorts = 0;

800
801
802
  if (clock) TIMER_TOC(timer_dosort);
}

803
/**
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
 * @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? */
822
  if (!cell_is_active_gravity(c, e)) return;
823
824

  /* Reset the gravity acceleration tensors */
825
  gravity_field_tensors_init(&c->multipole->pot, e->ti_current);
826
827
828
829
830
831
832
833
834
835
836

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

837
/**
838
839
840
841
842
 * @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.
843
 * @param timer Are we timing this ?
844
 */
845
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
846

847
#ifdef EXTRA_HYDRO_LOOP
848

849
  struct part *restrict parts = c->parts;
850
  struct xpart *restrict xparts = c->xparts;
851
  const int count = c->count;
852
  const struct engine *e = r->e;
853
  const struct cosmology *cosmo = e->cosmology;
854

855
856
  TIMER_TIC;

857
  /* Anything to do here? */
858
  if (!cell_is_active_hydro(c, e)) return;
859

860
861
862
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
863
      if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k], 0);
864
865
866
867
868
869
870
  } 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];
871
      struct xpart *restrict xp = &xparts[i];
872

873
      if (part_is_active(p, e)) {
874

875
        /* Finish the gradient calculation */
876
        hydro_end_gradient(p);
877
878
879
880
881
882
883
884
885
886
887

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

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

        /* The particle force values are now set.  Do _NOT_
           try to read any particle density variables! */

        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
888
889
890
      }
    }
  }
891

892
893
  if (timer) TIMER_TOC(timer_do_extra_ghost);

894
895
#else
  error("SWIFT was not compiled with the extra hydro loop activated.");
896
#endif
897
}
898

899
/**
900
901
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
902
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
903
 * @param r The runner thread.
904
 * @param c The cell.
905
 * @param timer Are we timing this ?
906
 */
907
void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
908

Matthieu Schaller's avatar
Matthieu Schaller committed
909
910
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
911
  const struct engine *e = r->e;
912
  const struct space *s = e->s;
913
914
  const struct hydro_space *hs = &s->hs;
  const struct cosmology *cosmo = e->cosmology;
lhausamm's avatar
lhausamm committed
915
  const struct chemistry_global_data *chemistry = e->chemistry;
916
  const float hydro_h_max = e->hydro_properties->h_max;
917
918
919
  const float eps = e->hydro_properties->h_tolerance;
  const float hydro_eta_dim =
      pow_dimension(e->hydro_properties->eta_neighbours);
920
  const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
921
  int redo = 0, count = 0;
922

923
924
  TIMER_TIC;

925
  /* Anything to do here? */
926
  if (!cell_is_active_hydro(c, e)) return;
927

928
929
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
930
    for (int k = 0; k < 8; k++)
931
932
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k], 0);
  } else {
933

934
    /* Init the list of active particles that have to be updated. */
935
    int *pid = NULL;
936
    if ((pid = (int *)malloc(sizeof(int) * c->count)) == NULL)
937
      error("Can't allocate memory for pid.");
938
939
940
941
942
    for (int k = 0; k < c->count; k++)
      if (part_is_active(&parts[k], e)) {
        pid[count] = k;
        ++count;
      }
943

944
945
946
    /* While there are particles that need to be updated... */
    for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
         num_reruns++) {
947

948
949
      /* Reset the redo-count. */
      redo = 0;
950

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

954
        /* Get a direct pointer on the part. */
955
956
        struct part *p = &parts[pid[i]];
        struct xpart *xp = &xparts[pid[i]];
957

958
#ifdef SWIFT_DEBUG_CHECKS
959
        /* Is this part within the timestep? */
960
961
962
        if (!part_is_active(p, e)) error("Ghost applied to inactive particle");
#endif

963
964
965
966
967
        /* 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;
968
        int has_no_neighbours = 0;
969

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

972
973
974
          /* Flag that there were no neighbours */
          has_no_neighbours = 1;

975
976
          /* Double h and try again */
          h_new = 2.f * h_old;