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

74
75
76
77
#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
78
#define TASK_LOOP_FEEDBACK 4
79

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

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

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

103
104
105
106
107
108
109
/* 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

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

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

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

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

137
  struct spart *restrict sparts = c->stars.parts;
138
139
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
140
141
142
143
144
145
  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;
146
  int redo = 0, scount = 0;
147

Loic Hausammann's avatar
Loic Hausammann committed
148
149
  double h_max = c->stars.h_max;

150
151
152
  TIMER_TIC;

  /* Anything to do here? */
153
  if (c->stars.count == 0) return;
Loic Hausammann's avatar
Loic Hausammann committed
154
  if (!cell_is_active_stars(c, e)) return;
155
156
157
158

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

162
163
164
        /* update h_max */
        if (c->progeny[k]->stars.h_max > h_max)
          h_max = c->progeny[k]->stars.h_max;
165
      }
166
167
168
169
  } else {

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

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

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

      /* Loop over the remaining active parts in this cell. */
198
      for (int i = 0; i < scount; i++) {
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? */
205
206
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
207
208
209
#endif

        /* Get some useful values */
210
        const float h_init = h_0[i];
211
212
213
        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);
214

215
216
217
        float h_new;
        int has_no_neighbours = 0;

218
        if (sp->density.wcount == 0.f) { /* No neighbours case */
219
220
221
222
223
224

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

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

226
227
228
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
229
          stars_end_density(sp, cosmo);
230
231

          /* Compute one step of the Newton-Raphson scheme */
232
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
233
          const float n_target = stars_eta_dim;
234
235
          const float f = n_sum - n_target;
          const float f_prime =
236
237
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
238

239
          /* Improve the bisection bounds */
240
241
242
243
          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);
244
245
246
247
248
249
250

#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

251
          /* Skip if h is already h_max and we don't have enough neighbours */
252
253
254
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
255

256
257
            stars_reset_feedback(sp);

258
259
260
261
262
            /* Ok, we are done with this particle */
            continue;
          }

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

264
265
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
266
267
268
269
270
271
272
273
274
275
276
277

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

278
279
280
281
282
283
284
285
286
#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);
287
288
289
290

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
291
292
293
294
295
296
        }

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

          /* Ok, correct then */
297
298
299
300
301
302
303
304
305
306
307
308
309
310

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

          /* If below the absolute maximum, try again */
313
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
314
315
316

            /* Flag for another round of fun */
            sid[redo] = sid[i];
317
318
319
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
320
321
322
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
323
            stars_init_spart(sp);
324
325
326
327

            /* Off we go ! */
            continue;

328
329
330
331
332
333
          } 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) {
334
335

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
336
            sp->h = stars_h_max;
337
338
339

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
340
              stars_spart_has_no_neighbours(sp, cosmo);
341
            }
342
343
344
345
346

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
347
348
349
          }
        }

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

352
353
        /* Check if h_max is increased */
        if (h_max < sp->h) h_max = sp->h;
Loic Hausammann's avatar
Loic Hausammann committed
354

355
        stars_reset_feedback(sp);
356

Loic Hausammann's avatar
Loic Hausammann committed
357
        /* Compute the stellar evolution  */
358
        stars_evolve_spart(sp, e->stars_properties, cosmo);
359
360
361
362
363
364
      }

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

      /* Re-set the counter for the next loop (potentially). */
365
366
      scount = redo;
      if (scount > 0) {
367
368
369
370
371

        /* 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
372
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
373
374
375
376
377
378
379
380

#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)
381
382
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
383
384
385
386
387
388

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

              /* Left or right? */
              if (l->t->ci == finger)
389
390
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
391
              else
392
393
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
394
395
396
397
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
398
399
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
400
401
402
403
404
405

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

              /* Left or right? */
              if (l->t->ci == finger)
406
407
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
408
              else
409
410
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
411
412
413
414
415
416
            }
          }
        }
      }
    }

417
418
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
419
420
421
    }

    /* Be clean */
422
423
    free(left);
    free(right);
424
    free(sid);
425
    free(h_0);
426
427
  }

428
  /* update h_max */
429
  if (h_max > c->stars.h_max) c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
430

431
432
433
  /* The ghost may not always be at the top level.
   * Therefore we need to update h_max above the ghost */
  if (c->stars.ghost) {
434
435
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
      if (h_max > tmp->stars.h_max) tmp->stars.h_max = h_max;
436
437
438
    }
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
439
  if (timer) TIMER_TOC(timer_dostars_ghost);
440
441
}

Tom Theuns's avatar
Tom Theuns committed
442
443
444
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
445
446
447
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
448
 */
449
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
450

451
452
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
453
454
455
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
456
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
457

458
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
459

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

Tom Theuns's avatar
Tom Theuns committed
463
464
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
465
    for (int k = 0; k < 8; k++)
466
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
467
  } else {
468

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

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

475
      /* Is this part within the time step? */
476
      if (gpart_is_active(gp, e)) {
477
478
        external_gravity_acceleration(time, potential, constants, gp);
      }
479
    }
480
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
481

482
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
483
484
}

485
486
487
488
489
490
491
492
493
/**
 * @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) {

494
495
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
  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
520
/**
521
522
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
523
524
525
526
527
528
529
 *
 * @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) {

530
  const struct engine *e = r->e;
531
532
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
533
534
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
535
  const struct unit_system *us = e->internal_units;
536
  const struct hydro_props *hydro_props = e->hydro_properties;
537
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
538
  const double time_base = e->time_base;
539
  const integertime_t ti_current = e->ti_current;
540
541
542
  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
543
544
545

  TIMER_TIC;

546
  /* Anything to do here? */
547
  if (!cell_is_active_hydro(c, e)) return;
548

Stefan Arridge's avatar
Stefan Arridge committed
549
550
551
552
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
553
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
554

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

558
559
560
      /* 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
561

562
      if (part_is_active(p, e)) {
563

564
        double dt_cool, dt_therm;
565
566
567
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
568
569
              get_integer_time_begin(ti_current - 1, p->time_bin);

570
571
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
572
573
574
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

575
576
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
577
          dt_therm = get_timestep(p->time_bin, time_base);
578
        }
579

580
        /* Let's cool ! */
581
582
583
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
584
      }
Stefan Arridge's avatar
Stefan Arridge committed
585
586
587
588
589
590
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
591
592
593
594
595
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

596
  struct engine *e = r->e;
597
  const struct cosmology *cosmo = e->cosmology;
598
599
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
600
601
602
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
603
  const int with_cosmology = (e->policy & engine_policy_cosmology);
604
  const int with_feedback = (e->policy & engine_policy_feedback);
605
606
607
  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;
608
  const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
609
610
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
611
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
612
613
614

  TIMER_TIC;

615
616
617
618
619
#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
620
  /* Anything to do here? */
621
  if (c->hydro.count == 0) return;
Matthieu Schaller's avatar
Matthieu Schaller committed
622
623
624
625
626
627
628
  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 {
629
630
631
632
633
634
635
636

    /* 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];

637
      /* Only work on active particles */
638
639
      if (part_is_active(p, e)) {

640
641
        /* 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
642
643
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
644

645
646
647
648
649
650
          /* 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);
651

652
653
654
655
656
657
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

          } else {
            dt_star = get_timestep(p->time_bin, time_base);
          }
658

659
660
661
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
662

663
664
665
666
667
668
669
          /* 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);

670
671
672
673
674
675
676
            /* 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);
            }
677
678
679
680
681
682
683
684
685
686
687
          }

        } 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
688
689
  }

690
691
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
692
  if (with_feedback && (c == c->top) &&
693
      (current_stars_count != c->stars.count)) {
694
695

    cell_clear_stars_sort_flags(c);
696
    runner_do_all_stars_sort(r, c);
697
  }
698

Matthieu Schaller's avatar
Matthieu Schaller committed
699
700
701
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
702
703
704
705
706
707
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
708
void runner_do_sort_ascending(struct entry *sort, int N) {
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762

  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
763
        }
764
765
766
767
768
769
770
771
772
773
774
775
      } 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
776
    }
777
778
779
  }
}

780
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
781
782
783
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
784
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
785
 */
786
787
788
789
790
791
792
793
794
#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
795
#else
Loic Hausammann's avatar
Loic Hausammann committed
796
797
798
799
#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
800
#endif
Loic Hausammann's avatar
Loic Hausammann committed
801
802
803

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
804

Pedro Gonnet's avatar
Pedro Gonnet committed
805
806
807
808
809
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
810
 * @param flags Cell flag.
811
812
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
813
814
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
815
 */
Loic Hausammann's avatar
Loic Hausammann committed
816
817
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
818
819

  struct entry *fingers[8];
820
821
822
  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
823
  float buff[8];
824

825
826
  TIMER_TIC;

827
828
829
830
#ifdef SWIFT_DEBUG_CHECKS
  if (c->hydro.super == NULL) error("Task called above the super level!!!");
#endif

831
  /* We need to do the local sorts plus whatever was requested further up. */
832
  flags |= c->hydro.do_sort;
833
  if (cleanup) {
834
    c->hydro.sorted = 0;
835
  } else {
836
    flags &= ~c->hydro.sorted;
837
  }
838
  if (flags == 0 && !c->hydro.do_sub_sort) return;
839
840

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

844
845
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
846
  runner_check_sorts_hydro(c, c->hydro.sorted);
847
848

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
849
850
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
851
852
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
853
  }
854

855
856
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
857
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
858
#endif
859

860
861
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
862
863
864
865
866

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

    /* Fill in the gaps within the progeny. */
867
    float dx_max_sort = 0.0f;
868
    float dx_max_sort_old = 0.0f;
869
    for (int k = 0; k < 8; k++) {
870
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
871
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
872
        runner_do_hydro_sort(r, c->progeny[k], flags,
873
874
875
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
876
877
878
        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);
879
      }
880
    }
881
882
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
883
884

    /* Loop over the 13 different sort arrays. */
885
    for (int j = 0; j < 13; j++) {
886
887
888
889
890

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

      /* Init the particle index offsets. */
891
      int off[8];
892
893
      off[0] = 0;
      for (int k = 1; k < 8; k++)
894
        if (c->progeny[k - 1] != NULL)
895
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
896
897
898
899
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
900
      int inds[8];
901
      for (int k = 0; k < 8; k++) {
902
        inds[k] = k;
903
904
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
905
906
907
908
909
910
911
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
912
913
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
914
          if (buff[inds[k]] < buff[inds[i]]) {
915
            int temp_i = inds[i];
916
917
918
919
920
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
921
      struct entry *finger = c->hydro.sort[j];
922
      for (int ind = 0; ind < count; ind++) {
923
924
925
926
927
928
929
930
931
932

        /* 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. */
933
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
934
          int temp_i = inds[k - 1];
935
936
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
937
        }
938

939
940
941
      } /* Merge. */

      /* Add a sentinel. */
942
943
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
944
945

      /* Mark as sorted. */
946
      atomic_or(&c->hydro.sorted, 1 << j);
947
948
949
950
951
952
953
954

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

955
    /* Reset the sort distance */
956
    if (c->hydro.sorted == 0) {
957
958
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
959
        error("Have non-NULL xparts in foreign cell");
960
#endif
961
962
963
964
965
966
967
968

      /* 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;
        }
969
      }
970
971
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
972
973
    }

974
    /* Fill the sort array. */
975
    for (int k = 0; k < count; k++) {
976
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
977
      for (int j = 0; j < 13; j++)