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
Loic Hausammann's avatar
Loic Hausammann committed
77
#define TASK_LOOP_FEEDBACK 4
78

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

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

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

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

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

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

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

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

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

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

149
150
151
  TIMER_TIC;

  /* Anything to do here? */
Loic Hausammann's avatar
Loic Hausammann committed
152
  if (!cell_is_active_stars(c, e)) return;
153
154
155
156

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

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

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

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

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

      /* Loop over the remaining active parts in this cell. */
196
      for (int i = 0; i < scount; i++) {
197
198
199
200
201
202

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

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

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

213
214
215
        float h_new;
        int has_no_neighbours = 0;

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

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

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

224
225
226
        } else {

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

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

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

#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

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

254
255
            stars_reset_feedback(sp);

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

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

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

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

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

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

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

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

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

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

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

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

            /* Off we go ! */
            continue;

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

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

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

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

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

	/* Check if h_max is increased */
	if (h_max < sp->h)
	  h_max = sp->h;

	stars_reset_feedback(sp);
355

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

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

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

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

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

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

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

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

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

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

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

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
431
  if (timer) TIMER_TOC(timer_dostars_ghost);
432
433
}

Tom Theuns's avatar
Tom Theuns committed
434
435
436
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
437
438
439
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
440
 */
441
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
442

443
444
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
445
446
447
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
448
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
449

450
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
451

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

Tom Theuns's avatar
Tom Theuns committed
455
456
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
457
    for (int k = 0; k < 8; k++)
458
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
459
  } else {
460

461
462
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
463

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

467
      /* Is this part within the time step? */
468
      if (gpart_is_active(gp, e)) {
469
470
        external_gravity_acceleration(time, potential, constants, gp);
      }
471
    }
472
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
473

474
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
475
476
}

477
478
479
480
481
482
483
484
485
/**
 * @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) {

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

522
  const struct engine *e = r->e;
523
524
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
525
526
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
527
  const struct unit_system *us = e->internal_units;
528
  const struct hydro_props *hydro_props = e->hydro_properties;
529
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
530
  const double time_base = e->time_base;
531
  const integertime_t ti_current = e->ti_current;
532
533
534
  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
535
536
537

  TIMER_TIC;

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

Stefan Arridge's avatar
Stefan Arridge committed
541
542
543
544
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
545
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
546

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

550
551
552
      /* 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
553

554
      if (part_is_active(p, e)) {
555

556
        double dt_cool, dt_therm;
557
558
559
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
560
561
              get_integer_time_begin(ti_current - 1, p->time_bin);

562
563
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
564
565
566
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

567
568
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
569
          dt_therm = get_timestep(p->time_bin, time_base);
570
        }
571

572
        /* Let's cool ! */
573
574
575
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
576
      }
Stefan Arridge's avatar
Stefan Arridge committed
577
578
579
580
581
582
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
583
584
585
586
587
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

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

  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 {
615
616
617
618
619
620
621
622

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

623
      /* Only work on active particles */
624
625
      if (part_is_active(p, e)) {

626
627
        /* 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
628
629
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
630

631
632
633
634
635
636
          /* 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);
637

638
639
640
641
642
643
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

645
646
647
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
648

649
650
651
652
653
654
655
656
657
658
659
660
661
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);

            /* 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
670
671
  }

672
673
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
674
675
  if (with_feedback && (c == c->hydro.super) &&
      (current_stars_count != c->stars.count)) {
676
677
    cell_clear_stars_sort_flags(c, /*is_super=*/1);
    runner_do_stars_sort(r, c, 0x1FFF, /*cleanup=*/0, /*timer=*/0);
678
  }
679

Matthieu Schaller's avatar
Matthieu Schaller committed
680
681
682
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
683
684
685
686
687
688
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
689
void runner_do_sort_ascending(struct entry *sort, int N) {
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
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743

  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
744
        }
745
746
747
748
749
750
751
752
753
754
755
756
      } 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
757
    }
758
759
760
  }
}

761
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
762
763
764
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
765
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
766
 */
767
768
769
770
771
772
773
774
775
#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
776
#else
Loic Hausammann's avatar
Loic Hausammann committed
777
778
779
780
#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
781
#endif
Loic Hausammann's avatar
Loic Hausammann committed
782
783
784

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
785

Pedro Gonnet's avatar
Pedro Gonnet committed
786
787
788
789
790
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
791
 * @param flags Cell flag.
792
793
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
794
795
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
796
 */
Loic Hausammann's avatar
Loic Hausammann committed
797
798
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
799
800

  struct entry *fingers[8];
801
802
803
  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
804
  float buff[8];
805

806
807
  TIMER_TIC;

808
  /* We need to do the local sorts plus whatever was requested further up. */
809
  flags |= c->hydro.do_sort;
810
  if (cleanup) {
811
    c->hydro.sorted = 0;
812
  } else {
813
    flags &= ~c->hydro.sorted;
814
  }
815
  if (flags == 0 && !c->hydro.do_sub_sort) return;
816
817

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

821
822
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
823
  runner_check_sorts_hydro(c, c->hydro.sorted);
824
825

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
826
827
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
828
829
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
830
  }
831

832
833
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
834
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
835
#endif
836

837
838
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
839
840
841
    if ((flags & (1 << j)) && c->hydro.sort[j] == NULL) {
      if ((c->hydro.sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                                     (count + 1))) == NULL)
842
843
        error("Failed to allocate sort memory.");
    }
844
845
846
847
848
849
  }

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

    /* Fill in the gaps within the progeny. */
850
    float dx_max_sort = 0.0f;
851
    float dx_max_sort_old = 0.0f;
852
    for (int k = 0; k < 8; k++) {
853
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
854
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
855
        runner_do_hydro_sort(r, c->progeny[k], flags,
856
857
858
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
859
860
861
        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);
862
      }
863
    }
864
865
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
866
867

    /* Loop over the 13 different sort arrays. */
868
    for (int j = 0; j < 13; j++) {
869
870
871
872
873

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

      /* Init the particle index offsets. */
874
      int off[8];
875
876
      off[0] = 0;
      for (int k = 1; k < 8; k++)
877
        if (c->progeny[k - 1] != NULL)
878
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
879
880
881
882
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
883
      int inds[8];
884
      for (int k = 0; k < 8; k++) {
885
        inds[k] = k;
886
887
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
888
889
890
891
892
893
894
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
895
896
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
897
          if (buff[inds[k]] < buff[inds[i]]) {
898
            int temp_i = inds[i];
899
900
901
902
903
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
904
      struct entry *finger = c->hydro.sort[j];
905
      for (int ind = 0; ind < count; ind++) {
906
907
908
909
910
911
912
913
914
915

        /* 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. */
916
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
917
          int temp_i = inds[k - 1];
918
919
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
920
        }
921

922
923
924
      } /* Merge. */

      /* Add a sentinel. */
925
926
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
927
928

      /* Mark as sorted. */
929
      atomic_or(&c->hydro.sorted, 1 << j);
930
931
932
933
934
935
936
937

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

938
    /* Reset the sort distance */
939
    if (c->hydro.sorted == 0) {
940
941
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
942
        error("Have non-NULL xparts in foreign cell");
943
#endif
944
945
946
947
948
949
950
951

      /* 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;
        }
952
      }
953
954
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
955
956
    }

957
    /* Fill the sort array. */
958
    for (int k = 0; k < count; k++) {
959
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
960
      for (int j = 0; j < 13; j++)
961
        if (flags & (1 << j)) {
962
963
964
965
          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];
966
        }
967
    }
968
969

    /* Add the sentinel and sort. */
970
    for (int j = 0; j < 13; j++)
971
      if (flags & (1 << j)) {
972
973
974
975
        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);
976
977
978
      }
  }

979
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
980
  /* Verify the sorting. */
981
  for (int j = 0; j < 13; j++) {
982
    if (!(flags & (1 << j))) continue;
983
    struct entry *finger = c->hydro.sort[j];
984
    for (int k = 1; k < count; k++) {
985
986
987
988
989