runner.c 109 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 "minmax.h"
James Willis's avatar
James Willis committed
59
#include "runner_doiact_vec.h"
60
#include "scheduler.h"
61
#include "sort_part.h"
62
#include "space.h"
63
#include "space_getsid.h"
64
#include "star_formation.h"
65
#include "stars.h"
66
67
#include "task.h"
#include "timers.h"
68
#include "timestep.h"
69
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
70
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
71

72
73
74
75
#define TASK_LOOP_DENSITY 0
#define TASK_LOOP_GRADIENT 1
#define TASK_LOOP_FORCE 2
#define TASK_LOOP_LIMITER 3
76

77
/* Import the density loop functions. */
78
#define FUNCTION density
79
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
80
#include "runner_doiact.h"
81
82
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
83

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

93
/* Import the force loop functions. */
94
#define FUNCTION force
95
#define FUNCTION_TASK_LOOP TASK_LOOP_FORCE
96
#include "runner_doiact.h"
97
98
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
99

100
101
102
103
104
105
106
/* 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

107
/* Import the gravity loop functions. */
108
#include "runner_doiact_grav.h"
109

110
111
/* Import the stars density loop functions. */
#define FUNCTION density
Loic Hausammann's avatar
Loic Hausammann committed
112
#define UPDATE_STARS 1
Loic Hausammann's avatar
Loic Hausammann committed
113
#include "runner_doiact_stars.h"
Loic Hausammann's avatar
Loic Hausammann committed
114
#undef UPDATE_STARS
115
116
117
118
#undef FUNCTION

/* Import the stars feedback loop functions. */
#define FUNCTION feedback
Loic Hausammann's avatar
Loic Hausammann committed
119
#include "runner_doiact_stars.h"
120
#undef FUNCTION
121

122
123
int emergency_sort = 0;

124
125
126
127
128
129
130
131
/**
 * @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
132
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
133

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

  TIMER_TIC;

  /* Anything to do here? */
Loic Hausammann's avatar
Loic Hausammann committed
148
  if (!cell_is_active_stars(c, e)) return;
149
150
151
152

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
Loic Hausammann's avatar
Loic Hausammann committed
153
      if (c->progeny[k] != NULL) runner_do_stars_ghost(r, c->progeny[k], 0);
154
155
156
157
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
158
159
160
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
161
    if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
Loic Hausammann's avatar
Loic Hausammann committed
162
      error("Can't allocate memory for sid.");
163
164
165
166
167
168
    if ((h_0 = (float *)malloc(sizeof(float) * c->hydro.count)) == NULL)
      error("Can't allocate memory for h_0.");
    if ((left = (float *)malloc(sizeof(float) * c->hydro.count)) == NULL)
      error("Can't allocate memory for left.");
    if ((right = (float *)malloc(sizeof(float) * c->hydro.count)) == NULL)
      error("Can't allocate memory for right.");
169
    for (int k = 0; k < c->stars.count; k++)
170
      if (spart_is_active(&sparts[k], e)) {
171
        sid[scount] = k;
172
173
174
        h_0[scount] = sparts[k].h;
        left[scount] = 0.f;
        right[scount] = stars_h_max;
175
        ++scount;
176
177
178
      }

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

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

      /* Loop over the remaining active parts in this cell. */
186
      for (int i = 0; i < scount; i++) {
187
188
189
190
191
192

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

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

        /* Get some useful values */
198
        const float h_init = h_0[i];
199
200
201
        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);
202

203
204
205
        float h_new;
        int has_no_neighbours = 0;

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

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

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

214
215
216
        } else {

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

          /* Compute one step of the Newton-Raphson scheme */
220
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
221
          const float n_target = stars_eta_dim;
222
223
          const float f = n_sum - n_target;
          const float f_prime =
224
225
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
226

227
228
229
230
231
232
233
234
235
236
          /* Improve the bisection bounds */
          if (n_sum < n_target) left[i] = max(left[i], h_old);
          if (n_sum > n_target) right[i] = min(right[i], h_old);

#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

237
          /* Skip if h is already h_max and we don't have enough neighbours */
238
239
240
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
241

242
243
            stars_reset_feedback(sp);

244
245
246
247
248
            /* Ok, we are done with this particle */
            continue;
          }

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

250
251
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
252
253
254
255
256
257
258
259
260
261
262
263

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

264
265
266
#ifdef SWIFT_DEBUG_CHECKS
          if ((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old))
            error(
267
                "Smoothing length correction not going in the right direction");
268
269
270
271
272
#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);
273
274
275
276

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
277
278
279
280
281
282
        }

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

          /* Ok, correct then */
283
284
285
286
287
288
289
290
291
292
293
294
295
296

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

          /* If below the absolute maximum, try again */
299
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
300
301
302

            /* Flag for another round of fun */
            sid[redo] = sid[i];
303
304
305
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
306
307
308
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
309
            stars_init_spart(sp);
310
311
312
313

            /* Off we go ! */
            continue;

314
315
316
317
318
319
          } 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) {
320
321

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
322
            sp->h = stars_h_max;
323
324
325

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
326
              stars_spart_has_no_neighbours(sp, cosmo);
327
            }
328
329
330
331
332

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
333
334
335
          }
        }

336
        /* We now have a particle whose smoothing length has converged */
337
        stars_reset_feedback(sp);
338

Loic Hausammann's avatar
Loic Hausammann committed
339
        /* Compute the stellar evolution  */
340
        stars_evolve_spart(sp, e->stars_properties, cosmo);
341
342
343
344
345
346
      }

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

      /* Re-set the counter for the next loop (potentially). */
347
348
      scount = redo;
      if (scount > 0) {
349
350
351
352
353

        /* 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
354
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
355
356
357
358
359
360
361
362

#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)
363
364
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
365
366
367
368
369
370

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

              /* Left or right? */
              if (l->t->ci == finger)
371
372
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
373
              else
374
375
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
376
377
378
379
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
380
381
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
382
383
384
385
386
387

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

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

399
400
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
401
402
403
    }

    /* Be clean */
404
405
    free(left);
    free(right);
406
    free(sid);
407
    free(h_0);
408
409
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
410
  if (timer) TIMER_TOC(timer_dostars_ghost);
411
412
}

Tom Theuns's avatar
Tom Theuns committed
413
414
415
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
416
417
418
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
419
 */
420
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
421

422
423
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
424
425
426
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
427
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
428

429
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
430

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

Tom Theuns's avatar
Tom Theuns committed
434
435
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
436
    for (int k = 0; k < 8; k++)
437
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
438
  } else {
439

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

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

446
      /* Is this part within the time step? */
447
      if (gpart_is_active(gp, e)) {
448
449
        external_gravity_acceleration(time, potential, constants, gp);
      }
450
    }
451
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
452

453
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
454
455
}

456
457
458
459
460
461
462
463
464
/**
 * @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) {

465
466
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
  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
491
/**
492
493
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
494
495
496
497
498
499
500
 *
 * @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) {

501
  const struct engine *e = r->e;
502
503
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
504
505
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
506
  const struct unit_system *us = e->internal_units;
507
  const struct hydro_props *hydro_props = e->hydro_properties;
508
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
509
  const double time_base = e->time_base;
510
  const integertime_t ti_current = e->ti_current;
511
512
513
  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
514
515
516

  TIMER_TIC;

517
  /* Anything to do here? */
518
  if (!cell_is_active_hydro(c, e)) return;
519

Stefan Arridge's avatar
Stefan Arridge committed
520
521
522
523
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
524
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
525

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

529
530
531
      /* 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
532

533
      if (part_is_active(p, e)) {
534

535
        double dt_cool, dt_therm;
536
537
538
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
539
540
              get_integer_time_begin(ti_current - 1, p->time_bin);

541
542
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
543
544
545
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

546
547
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
548
          dt_therm = get_timestep(p->time_bin, time_base);
549
        }
550

551
        /* Let's cool ! */
552
553
554
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
555
      }
Stefan Arridge's avatar
Stefan Arridge committed
556
557
558
559
560
561
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
562
563
564
/**
 *
 */
565
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
Matthieu Schaller's avatar
Matthieu Schaller committed
566

567
  struct engine *e = r->e;
568
  const struct cosmology *cosmo = e->cosmology;
569
570
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
571
572
573
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
574
  const int with_cosmology = (e->policy & engine_policy_cosmology);
575
576
577
  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;
578
579
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
580
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
581
582
583
584
585
586
587
588
589

  TIMER_TIC;

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

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

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

600
      /* Only work on active particles */
601
602
      if (part_is_active(p, e)) {

603
604
605
        /* Is this particle star forming? */
        if (star_formation_is_star_forming(p, xp, sf_props, phys_const, cosmo,
                                           hydro_props, us, cooling)) {
606

607
608
609
610
611
612
          /* 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);
613

614
615
616
617
618
619
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

621
622
623
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
624

625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
          /* 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);

            /* Copy the properties of the gas particle to the star particle */
            star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
                                           with_cosmology);
          }

        } 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
646
647
  }

648
649
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
650
651
652
  if ((c == c->hydro.super) && (current_stars_count != c->stars.count)) {
    cell_clear_stars_sort_flags(c, /*is_super=*/1);
    runner_do_stars_sort(r, c, 0x1FFF, /*cleanup=*/0, /*timer=*/0);
653
  }
654

Matthieu Schaller's avatar
Matthieu Schaller committed
655
656
657
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
658
659
660
661
662
663
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
664
void runner_do_sort_ascending(struct entry *sort, int N) {
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718

  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
719
        }
720
721
722
723
724
725
726
727
728
729
730
731
      } 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
732
    }
733
734
735
  }
}

736
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
737
738
739
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
740
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
741
 */
742
743
744
745
746
747
748
749
750
#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
751
#else
Loic Hausammann's avatar
Loic Hausammann committed
752
753
754
755
#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
756
#endif
Loic Hausammann's avatar
Loic Hausammann committed
757
758
759

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
760

Pedro Gonnet's avatar
Pedro Gonnet committed
761
762
763
764
765
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
766
 * @param flags Cell flag.
767
768
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
769
770
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
771
 */
Loic Hausammann's avatar
Loic Hausammann committed
772
773
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
774
775

  struct entry *fingers[8];
776
777
778
  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
779
  float buff[8];
780

781
782
  TIMER_TIC;

783
  /* We need to do the local sorts plus whatever was requested further up. */
784
  flags |= c->hydro.do_sort;
785
  if (cleanup) {
786
    c->hydro.sorted = 0;
787
  } else {
788
    flags &= ~c->hydro.sorted;
789
  }
790
  if (flags == 0 && !c->hydro.do_sub_sort) return;
791
792

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

796
797
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
798
  runner_check_sorts_hydro(c, c->hydro.sorted);
799
800

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
801
802
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
803
804
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
805
  }
806

807
808
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
809
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
810
#endif
811

812
813
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
814
815
816
    if ((flags & (1 << j)) && c->hydro.sort[j] == NULL) {
      if ((c->hydro.sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                                     (count + 1))) == NULL)
817
818
        error("Failed to allocate sort memory.");
    }
819
820
821
822
823
824
  }

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

    /* Fill in the gaps within the progeny. */
825
    float dx_max_sort = 0.0f;
826
    float dx_max_sort_old = 0.0f;
827
    for (int k = 0; k < 8; k++) {
828
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
829
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
830
        runner_do_hydro_sort(r, c->progeny[k], flags,
831
832
833
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
834
835
836
        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);
837
      }
838
    }
839
840
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
841
842

    /* Loop over the 13 different sort arrays. */
843
    for (int j = 0; j < 13; j++) {
844
845
846
847
848

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

      /* Init the particle index offsets. */
849
      int off[8];
850
851
      off[0] = 0;
      for (int k = 1; k < 8; k++)
852
        if (c->progeny[k - 1] != NULL)
853
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
854
855
856
857
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
858
      int inds[8];
859
      for (int k = 0; k < 8; k++) {
860
        inds[k] = k;
861
862
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
863
864
865
866
867
868
869
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
870
871
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
872
          if (buff[inds[k]] < buff[inds[i]]) {
873
            int temp_i = inds[i];
874
875
876
877
878
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
879
      struct entry *finger = c->hydro.sort[j];
880
      for (int ind = 0; ind < count; ind++) {
881
882
883
884
885
886
887
888
889
890

        /* 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. */
891
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
892
          int temp_i = inds[k - 1];
893
894
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
895
        }
896

897
898
899
      } /* Merge. */

      /* Add a sentinel. */
900
901
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
902
903

      /* Mark as sorted. */
904
      atomic_or(&c->hydro.sorted, 1 << j);
905
906
907
908
909
910
911
912

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

913
    /* Reset the sort distance */
914
    if (c->hydro.sorted == 0) {
915
916
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
917
        error("Have non-NULL xparts in foreign cell");
918
#endif
919
920
921
922
923
924
925
926

      /* 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;
        }
927
      }
928
929
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
930
931
    }

932
    /* Fill the sort array. */
933
    for (int k = 0; k < count; k++) {
934
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
935
      for (int j = 0; j < 13; j++)
936
        if (flags & (1 << j)) {
937
938
939
940
          c->hydro.sort[j][k].i = k;
          c->hydro.sort[j][k].d = px[0] * runner_shift[j][0] +
                                  px[1] * runner_shift[j][1] +
                                  px[2] * runner_shift[j][2];
941
        }
942
    }
943
944

    /* Add the sentinel and sort. */
945
    for (int j = 0; j < 13; j++)
946
      if (flags & (1 << j)) {
947
948
949
950
        c->hydro.sort[j][count].d = FLT_MAX;
        c->hydro.sort[j][count].i = 0;
        runner_do_sort_ascending(c->hydro.sort[j], count);
        atomic_or(&c->hydro.sorted, 1 << j);
951
952
953
      }
  }

954
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
955
  /* Verify the sorting. */
956
  for (int j = 0; j < 13; j++) {
957
    if (!(flags & (1 << j))) continue;
958
    struct entry *finger = c->hydro.sort[j];
959
    for (int k = 1; k < count; k++) {
960
961
962
963
964
      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
965

966
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
967
  runner_check_sorts_hydro(c, flags);
968
969

  /* Make sure the sort flags are consistent (upward). */
Pedro Gonnet's avatar
Pedro Gonnet committed
970
971
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
972
973
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags.");
974
  }
975
#endif
976

977
  /* Clear the cell's sort flags. */
978
979
980
  c->hydro.do_sort = 0;
  c->hydro.do_sub_sort = 0;
  c->hydro.requires_sorts = 0;
981