runner.c 116 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 "entropy_floor.h"
52
#include "error.h"
53
#include "feedback.h"
54
55
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
56
#include "hydro_properties.h"
57
#include "kick.h"
Loikki's avatar
Loikki committed
58
#include "logger.h"
59
#include "memuse.h"
60
#include "minmax.h"
James Willis's avatar
James Willis committed
61
#include "runner_doiact_vec.h"
62
#include "scheduler.h"
63
#include "sort_part.h"
64
#include "space.h"
65
#include "space_getsid.h"
66
#include "star_formation.h"
67
#include "star_formation_iact.h"
68
#include "stars.h"
69
70
#include "task.h"
#include "timers.h"
71
#include "timestep.h"
72
#include "timestep_limiter.h"
73
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
74

75
76
77
78
#define TASK_LOOP_DENSITY 0
#define TASK_LOOP_GRADIENT 1
#define TASK_LOOP_FORCE 2
#define TASK_LOOP_LIMITER 3
Loic Hausammann's avatar
Loic Hausammann committed
79
#define TASK_LOOP_FEEDBACK 4
80

81
/* Import the density loop functions. */
82
#define FUNCTION density
83
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
84
#include "runner_doiact.h"
85
86
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
87

88
/* Import the gradient loop functions (if required). */
89
90
#ifdef EXTRA_HYDRO_LOOP
#define FUNCTION gradient
91
#define FUNCTION_TASK_LOOP TASK_LOOP_GRADIENT
92
#include "runner_doiact.h"
93
94
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
95
96
#endif

97
/* Import the force loop functions. */
98
#define FUNCTION force
99
#define FUNCTION_TASK_LOOP TASK_LOOP_FORCE
100
#include "runner_doiact.h"
101
102
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
103

104
105
106
107
108
109
110
/* Import the limiter loop functions. */
#define FUNCTION limiter
#define FUNCTION_TASK_LOOP TASK_LOOP_LIMITER
#include "runner_doiact.h"
#undef FUNCTION
#undef FUNCTION_TASK_LOOP

111
/* Import the gravity loop functions. */
112
#include "runner_doiact_grav.h"
113

114
115
/* Import the stars density loop functions. */
#define FUNCTION density
Loic Hausammann's avatar
Loic Hausammann committed
116
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
Loic Hausammann's avatar
Loic Hausammann committed
117
#include "runner_doiact_stars.h"
Loic Hausammann's avatar
Loic Hausammann committed
118
#undef FUNCTION_TASK_LOOP
119
120
121
122
#undef FUNCTION

/* Import the stars feedback loop functions. */
#define FUNCTION feedback
Loic Hausammann's avatar
Loic Hausammann committed
123
#define FUNCTION_TASK_LOOP TASK_LOOP_FEEDBACK
Loic Hausammann's avatar
Loic Hausammann committed
124
#include "runner_doiact_stars.h"
Loic Hausammann's avatar
Loic Hausammann committed
125
#undef FUNCTION_TASK_LOOP
126
#undef FUNCTION
127
128
129
130
131
132
133
134
135

/**
 * @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
136
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
137

138
  struct spart *restrict sparts = c->stars.parts;
139
  const struct engine *e = r->e;
140
141
  const struct unit_system *us = e->internal_units;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
142
  const struct cosmology *cosmo = e->cosmology;
143
  const struct feedback_props *feedback_props = e->feedback_props;
144
145
146
147
148
149
  const float stars_h_max = e->hydro_properties->h_max;
  const float stars_h_min = e->hydro_properties->h_min;
  const float eps = e->stars_properties->h_tolerance;
  const float stars_eta_dim =
      pow_dimension(e->stars_properties->eta_neighbours);
  const int max_smoothing_iter = e->stars_properties->max_smoothing_iterations;
150
  int redo = 0, scount = 0;
151

Matthieu Schaller's avatar
Matthieu Schaller committed
152
  /* Running value of the maximal smoothing length */
Loic Hausammann's avatar
Loic Hausammann committed
153
154
  double h_max = c->stars.h_max;

155
156
157
  TIMER_TIC;

  /* Anything to do here? */
158
  if (c->stars.count == 0) return;
Loic Hausammann's avatar
Loic Hausammann committed
159
  if (!cell_is_active_stars(c, e)) return;
160
161
162

  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
163
    for (int k = 0; k < 8; k++) {
164
      if (c->progeny[k] != NULL) {
165
        runner_do_stars_ghost(r, c->progeny[k], 0);
166

Matthieu Schaller's avatar
Matthieu Schaller committed
167
168
        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->stars.h_max);
169
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
170
    }
171
172
173
174
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
175
176
177
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
178
    if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
Loic Hausammann's avatar
Loic Hausammann committed
179
      error("Can't allocate memory for sid.");
180
    if ((h_0 = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
181
      error("Can't allocate memory for h_0.");
182
    if ((left = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
183
      error("Can't allocate memory for left.");
184
    if ((right = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
185
      error("Can't allocate memory for right.");
186
    for (int k = 0; k < c->stars.count; k++)
187
      if (spart_is_active(&sparts[k], e)) {
188
        sid[scount] = k;
189
190
191
        h_0[scount] = sparts[k].h;
        left[scount] = 0.f;
        right[scount] = stars_h_max;
192
        ++scount;
193
194
195
      }

    /* While there are particles that need to be updated... */
196
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
197
         num_reruns++) {
198

199
200
201
202
      /* Reset the redo-count. */
      redo = 0;

      /* Loop over the remaining active parts in this cell. */
203
      for (int i = 0; i < scount; i++) {
204
205
206
207
208
209

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

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

        /* Get some useful values */
215
        const float h_init = h_0[i];
216
217
218
        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);
219

220
221
222
        float h_new;
        int has_no_neighbours = 0;

223
        if (sp->density.wcount == 0.f) { /* No neighbours case */
224
225
226
227
228
229

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

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

231
232
233
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
234
          stars_end_density(sp, cosmo);
235
236

          /* Compute one step of the Newton-Raphson scheme */
237
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
238
          const float n_target = stars_eta_dim;
239
240
          const float f = n_sum - n_target;
          const float f_prime =
241
242
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
243

244
          /* Improve the bisection bounds */
245
246
247
248
          if (n_sum < n_target)
            left[i] = max(left[i], h_old);
          else if (n_sum > n_target)
            right[i] = min(right[i], h_old);
249
250
251
252
253
254
255

#ifdef SWIFT_DEBUG_CHECKS
          /* Check the validity of the left and right bounds */
          if (left[i] > right[i])
            error("Invalid left (%e) and right (%e)", left[i], right[i]);
#endif

256
          /* Skip if h is already h_max and we don't have enough neighbours */
257
258
259
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
260

261
            stars_reset_feedback(sp);
262
263
264
265
266
267

            /* Ok, we are done with this particle */
            continue;
          }

          /* Normal case: Use Newton-Raphson to get a better value of h */
268

269
270
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
271
272
273
274
275
276
277
278
279
280
281
282

          /* Be verbose about the particles that struggle to converge */
          if (num_reruns > max_smoothing_iter - 10) {

            message(
                "Smoothing length convergence problem: iter=%d p->id=%lld "
                "h_init=%12.8e h_old=%12.8e h_new=%12.8e f=%f f_prime=%f "
                "n_sum=%12.8e n_target=%12.8e left=%12.8e right=%12.8e",
                num_reruns, sp->id, h_init, h_old, h_new, f, f_prime, n_sum,
                n_target, left[i], right[i]);
          }

283
284
285
286
287
288
289
290
291
#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);
292
293
294
295

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
296
297
298
299
300
301
        }

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

          /* Ok, correct then */
302
303
304
305
306
307
308
309
310
311
312
313
314
315

          /* Case where we have been oscillating around the solution */
          if ((h_new == left[i] && h_old == right[i]) ||
              (h_old == left[i] && h_new == right[i])) {

            /* Bissect the remaining interval */
            sp->h = pow_inv_dimension(
                0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));

          } else {

            /* Normal case */
            sp->h = h_new;
          }
316
317

          /* If below the absolute maximum, try again */
318
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
319
320
321

            /* Flag for another round of fun */
            sid[redo] = sid[i];
322
323
324
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
325
326
327
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
328
            stars_init_spart(sp);
329
330
331
332

            /* Off we go ! */
            continue;

333
334
335
336
337
338
          } else if (sp->h <= stars_h_min) {

            /* Ok, this particle is a lost cause... */
            sp->h = stars_h_min;

          } else if (sp->h >= stars_h_max) {
339
340

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
341
            sp->h = stars_h_max;
342
343
344

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
345
              stars_spart_has_no_neighbours(sp, cosmo);
346
            }
347
348
349
350
351

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
352
353
354
          }
        }

355
        /* We now have a particle whose smoothing length has converged */
Loic Hausammann's avatar
Loic Hausammann committed
356

Matthieu Schaller's avatar
Matthieu Schaller committed
357
358
        /* Check if h_max has increased */
        h_max = max(h_max, sp->h);
Loic Hausammann's avatar
Loic Hausammann committed
359

360
        stars_reset_feedback(sp);
361

362
        /* Only do feedback if stars have a reasonable birth time */
363
364
        if (sp->birth_time != -1.) {

365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
          const integertime_t ti_step = get_integer_timestep(sp->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(e->ti_current - 1, sp->time_bin);

          /* Get particle time-step */
          double dt;
          if (with_cosmology) {
            dt = cosmology_get_delta_time(e->cosmology, ti_begin,
                                          ti_begin + ti_step);
          } else {
            dt = get_timestep(sp->time_bin, e->time_base);
          }

          /* Calculate age of the star at current time */
          double star_age_end_of_step;
380
          if (with_cosmology) {
381
            star_age_end_of_step = cosmology_get_delta_time_from_scale_factors(
382
383
                cosmo, sp->birth_scale_factor, (float)cosmo->a);
          } else {
384
            star_age_end_of_step = e->time - sp->birth_time;
385
386
          }

387
388
389
390
391
392
393
394
395
396
397
398
399
400
          /* Reset the feedback fields of the star particle */
          feedback_prepare_spart(sp, feedback_props);

          /* Has this star been around for a while ? */
          if (star_age_end_of_step > 0.) {

            /* Age of the star at the start of the step */
            const double star_age_beg_of_step =
                max(star_age_end_of_step - dt, 0.);

            /* Compute the stellar evolution  */
            feedback_evolve_spart(sp, feedback_props, cosmo, us,
                                  star_age_beg_of_step, dt);
          }
401
        }
402
403
404
405
406
407
      }

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

      /* Re-set the counter for the next loop (potentially). */
408
409
      scount = redo;
      if (scount > 0) {
410
411
412
413
414

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

          /* Run through this cell's density interactions. */
Loic Hausammann's avatar
Loic Hausammann committed
415
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
416
417
418
419
420
421
422
423

#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)
424
425
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
426
427
428
429
430
431

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

              /* Left or right? */
              if (l->t->ci == finger)
432
433
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
434
              else
435
436
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
437
438
439
440
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
441
442
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
443
444
445
446
447
448

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

              /* Left or right? */
              if (l->t->ci == finger)
449
450
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
451
              else
452
453
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
454
455
456
457
458
459
            }
          }
        }
      }
    }

460
461
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
462
463
464
    }

    /* Be clean */
465
466
    free(left);
    free(right);
467
    free(sid);
468
    free(h_0);
469
470
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
471
472
  /* Update h_max */
  c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
473

474
  /* The ghost may not always be at the top level.
475
   * Therefore we need to update h_max between the super- and top-levels */
476
  if (c->stars.ghost) {
477
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
478
      atomic_max_d(&tmp->stars.h_max, h_max);
479
480
481
    }
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
482
  if (timer) TIMER_TOC(timer_dostars_ghost);
483
484
}

Tom Theuns's avatar
Tom Theuns committed
485
486
487
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
488
489
490
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
491
 */
492
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
493

494
495
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
496
497
498
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
499
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
500

501
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
502

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

Tom Theuns's avatar
Tom Theuns committed
506
507
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
508
    for (int k = 0; k < 8; k++)
509
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
510
  } else {
511

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

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

518
      /* Is this part within the time step? */
519
      if (gpart_is_active(gp, e)) {
520
521
        external_gravity_acceleration(time, potential, constants, gp);
      }
522
    }
523
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
524

525
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
526
527
}

528
529
530
531
532
533
534
535
536
/**
 * @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) {

537
538
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
  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
563
/**
564
565
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
566
567
568
569
570
571
572
 *
 * @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) {

573
  const struct engine *e = r->e;
574
575
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
576
577
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
578
  const struct unit_system *us = e->internal_units;
579
  const struct hydro_props *hydro_props = e->hydro_properties;
580
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
581
  const double time_base = e->time_base;
582
  const integertime_t ti_current = e->ti_current;
583
584
585
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
  const int count = c->hydro.count;
Stefan Arridge's avatar
Stefan Arridge committed
586
587
588

  TIMER_TIC;

589
  /* Anything to do here? */
590
  if (!cell_is_active_hydro(c, e)) return;
591

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

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

601
602
603
      /* 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
604

605
      if (part_is_active(p, e)) {
606

607
        double dt_cool, dt_therm;
608
609
610
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
611
612
              get_integer_time_begin(ti_current - 1, p->time_bin);

613
614
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
615
616
617
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

618
619
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
620
          dt_therm = get_timestep(p->time_bin, time_base);
621
        }
622

623
        /* Let's cool ! */
624
625
626
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
627
      }
Stefan Arridge's avatar
Stefan Arridge committed
628
629
630
631
632
633
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
634
635
636
637
638
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

639
  struct engine *e = r->e;
640
  const struct cosmology *cosmo = e->cosmology;
641
642
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
643
644
645
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
646
  const int with_cosmology = (e->policy & engine_policy_cosmology);
647
  const int with_feedback = (e->policy & engine_policy_feedback);
648
649
650
  const struct hydro_props *restrict hydro_props = e->hydro_properties;
  const struct unit_system *restrict us = e->internal_units;
  struct cooling_function_data *restrict cooling = e->cooling_func;
651
  const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
652
653
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
654
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
655
656
657

  TIMER_TIC;

658
659
660
661
662
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID != e->nodeID)
    error("Running star formation task on a foreign node!");
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
663
  /* Anything to do here? */
664
  if (c->hydro.count == 0) return;
Matthieu Schaller's avatar
Matthieu Schaller committed
665
666
667
668
669
670
671
  if (!cell_is_active_hydro(c, e)) return;

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_star_formation(r, c->progeny[k], 0);
  } else {
672
673
674
675
676
677
678
679

    /* Loop over the gas particles 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];

680
      /* Only work on active particles */
681
682
      if (part_is_active(p, e)) {

683
684
        /* Is this particle star forming? */
        if (star_formation_is_star_forming(p, xp, sf_props, phys_const, cosmo,
Folkert Nobels's avatar
Folkert Nobels committed
685
686
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
687

688
689
690
691
692
693
          /* Time-step size for this particle */
          double dt_star;
          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);
694

695
696
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
697

698
699
700
          } else {
            dt_star = get_timestep(p->time_bin, time_base);
          }
701

702
703
704
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
705

706
707
708
709
710
711
712
          /* Are we forming a star particle from this SF rate? */
          if (star_formation_should_convert_to_star(p, xp, sf_props, e,
                                                    dt_star)) {

            /* Convert the gas particle to a star particle */
            struct spart *sp = cell_convert_part_to_spart(e, c, p, xp);

713
714
715
716
717
718
719
            /* Did we get a star? (Or did we run out of spare ones?) */
            if (sp != NULL) {

              /* Copy the properties of the gas particle to the star particle */
              star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
                                             with_cosmology);
            }
720
721
722
723
724
725
726
727
728
729
730
          }

        } else { /* Are we not star-forming? */

          /* Update the particle to flag it as not star-forming */
          star_formation_update_part_not_SFR(p, xp, e, sf_props,
                                             with_cosmology);

        } /* Not Star-forming? */
      }   /* is active? */
    }     /* Loop over particles */
Matthieu Schaller's avatar
Matthieu Schaller committed
731
732
  }

733
734
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
735
  if (with_feedback && (c == c->top) &&
736
      (current_stars_count != c->stars.count)) {
737
738

    cell_clear_stars_sort_flags(c);
739
    runner_do_all_stars_sort(r, c);
740
  }
741

Matthieu Schaller's avatar
Matthieu Schaller committed
742
743
744
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
745
746
747
748
749
750
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
751
void runner_do_sort_ascending(struct entry *sort, int N) {
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805

  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
806
        }
807
808
809
810
811
812
813
814
815
816
817
818
      } 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
819
    }
820
821
822
  }
}

823
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
824
825
826
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
827
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
828
 */
829
830
831
832
833
834
835
836
837
#define RUNNER_CHECK_SORTS(TYPE)                                               \
  void runner_check_sorts_##TYPE(struct cell *c, int flags) {                  \
                                                                               \
    if (flags & ~c->TYPE.sorted) error("Inconsistent sort flags (downward)!"); \
    if (c->split)                                                              \
      for (int k = 0; k < 8; k++)                                              \
        if (c->progeny[k] != NULL && c->progeny[k]->TYPE.count > 0)            \
          runner_check_sorts_##TYPE(c->progeny[k], c->TYPE.sorted);            \
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
838
#else
Loic Hausammann's avatar
Loic Hausammann committed
839
840
841
842
#define RUNNER_CHECK_SORTS(TYPE)                                       \
  void runner_check_sorts_##TYPE(struct cell *c, int flags) {          \
    error("Calling debugging code without debugging flag activated."); \
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
843
#endif
Loic Hausammann's avatar
Loic Hausammann committed
844
845
846

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
847

Pedro Gonnet's avatar
Pedro Gonnet committed
848
849
850
851
852
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
853
 * @param flags Cell flag.
854
855
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
856
857
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
858
 */
Loic Hausammann's avatar
Loic Hausammann committed
859
860
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
861
862

  struct entry *fingers[8];
863
864
865
  const int count = c->hydro.count;
  const struct part *parts = c->hydro.parts;
  struct xpart *xparts = c->hydro.xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
866
  float buff[8];
867

868
869
  TIMER_TIC;

870
871
872
873
#ifdef SWIFT_DEBUG_CHECKS
  if (c->hydro.super == NULL) error("Task called above the super level!!!");
#endif

874
  /* We need to do the local sorts plus whatever was requested further up. */
875
  flags |= c->hydro.do_sort;
876
  if (cleanup) {
877
    c->hydro.sorted = 0;
878
  } else {
879
    flags &= ~c->hydro.sorted;
880
  }
881
  if (flags == 0 && !c->hydro.do_sub_sort) return;
882
883

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

887
888
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
889
  runner_check_sorts_hydro(c, c->hydro.sorted);
890
891

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
892
893
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
894
895
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
896
  }
897

898
899
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
900
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
901
#endif
902

903
904
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
905
906
907
908
909

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

    /* Fill in the gaps within the progeny. */
910
    float dx_max_sort = 0.0f;
911
    float dx_max_sort_old = 0.0f;
912
    for (int k = 0; k < 8; k++) {
913
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
914
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
915
        runner_do_hydro_sort(r, c->progeny[k], flags,
916
917
918
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
919
920
921
        dx_max_sort = max(dx_max_sort, c->progeny[k]->hydro.dx_max_sort);
        dx_max_sort_old =
            max(dx_max_sort_old, c->progeny[k]->hydro.dx_max_sort_old);
922
      }
923
    }
924
925
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
926
927

    /* Loop over the 13 different sort arrays. */
928
    for (int j = 0; j < 13; j++) {
929
930
931
932
933

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

      /* Init the particle index offsets. */
934
      int off[8];
935
936
      off[0] = 0;
      for (int k = 1; k < 8; k++)
937
        if (c->progeny[k - 1] != NULL)
938
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
939
940
941
942
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
943
      int inds[8];
944
      for (int k = 0; k < 8; k++) {
945
        inds[k] = k;
946
947
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
948
949
950
951
952
953
954
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
955
956
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
957
          if (buff[inds[k]] < buff[inds[i]]) {
958
            int temp_i = inds[i];
959
960
961
962
963
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
964
      struct entry *finger = c->hydro.sort[j];
965
      for (int ind = 0; ind < count; ind++) {
966
967
968
969
970
971
972
973
974
975

        /* 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. */
976
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
977
          int temp_i = inds[k - 1];
978
979
          inds[k - 1] = inds[k];
          inds[k] = temp_i;