runner.c 110 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 "star_formation_iact.h"
66
#include "stars.h"
67
68
#include "task.h"
#include "timers.h"
69
#include "timestep.h"
70
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
71
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
72

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

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

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

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

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

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

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

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

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

202
203
204
        float h_new;
        int has_no_neighbours = 0;

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

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

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

213
214
215
        } else {

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

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

226
          /* Improve the bisection bounds */
227
228
229
230
          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);
231
232
233
234
235
236
237

#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

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

243
244
            stars_reset_feedback(sp);

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

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

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

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

265
266
267
268
269
270
271
272
273
#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);
274
275
276
277

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

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

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

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

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

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

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

            /* Off we go ! */
            continue;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

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

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

  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++)
      if (c->progeny[k] != NULL) runner_do_star_formation(r, c->progeny[k], 0);
  } else {
595
596
597
598
599
600
601
602

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

603
      /* Only work on active particles */
604
605
      if (part_is_active(p, e)) {

606
607
        /* Is this particle star forming? */
        if (star_formation_is_star_forming(p, xp, sf_props, phys_const, cosmo,
608
                                           hydro_props, us, cooling, entropy)) {
609

610
611
612
613
614
615
          /* 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);
616

617
618
619
620
621
622
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

624
625
626
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
627

628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
          /* 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
649
650
  }

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

Matthieu Schaller's avatar
Matthieu Schaller committed
659
660
661
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
662
663
664
665
666
667
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
668
void runner_do_sort_ascending(struct entry *sort, int N) {
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
719
720
721
722

  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
723
        }
724
725
726
727
728
729
730
731
732
733
734
735
      } 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
736
    }
737
738
739
  }
}

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

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
764

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

  struct entry *fingers[8];
780
781
782
  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
783
  float buff[8];
784

785
786
  TIMER_TIC;

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

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

800
801
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
802
  runner_check_sorts_hydro(c, c->hydro.sorted);
803
804

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

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

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

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

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

    /* Loop over the 13 different sort arrays. */
847
    for (int j = 0; j < 13; j++) {
848
849
850
851
852

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

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

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

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

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

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

901
902
903
      } /* Merge. */

      /* Add a sentinel. */
904
905
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
906
907

      /* Mark as sorted. */
908
      atomic_or(&c->hydro.sorted, 1 << j);
909
910
911
912
913
914
915
916

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

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

      /* 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;
        }
931
      }
932
933
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
934
935
    }

936
    /* Fill the sort array. */
937
    for (int k = 0; k < count; k++) {
938
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
939
      for (int j = 0; j < 13; j++)
940
        if (flags & (1 << j)) {
941
942
943
944
          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];
945
        }
946
    }
947
948

    /* Add the sentinel and sort. */
949
    for (int j = 0; j < 13; j++)
950
      if (flags & (1 << j)) {
951
952
953
954
        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);
955
956
957
      }
  }

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

970
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
971
  runner_check_sorts_hydro(c, flags);
972
973

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