runner.c 122 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 "black_holes.h"
45
#include "cell.h"
46
#include "chemistry.h"
47
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
48
#include "cooling.h"
49
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
50
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
51
#include "engine.h"
52
#include "entropy_floor.h"
53
#include "error.h"
54
#include "feedback.h"
55
56
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
57
#include "hydro_properties.h"
58
#include "kick.h"
Loikki's avatar
Loikki committed
59
#include "logger.h"
60
#include "memuse.h"
61
#include "minmax.h"
James Willis's avatar
James Willis committed
62
#include "runner_doiact_vec.h"
63
#include "scheduler.h"
64
#include "sort_part.h"
65
#include "space.h"
66
#include "space_getsid.h"
67
#include "star_formation.h"
68
#include "star_formation_iact.h"
69
#include "star_formation_logger.h"
70
#include "stars.h"
71
72
#include "task.h"
#include "timers.h"
73
#include "timestep.h"
74
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
75
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
76

77
78
79
80
#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
81
#define TASK_LOOP_FEEDBACK 4
82

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

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

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

106
107
108
109
110
111
112
/* 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

113
/* Import the gravity loop functions. */
114
#include "runner_doiact_grav.h"
115

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

/* Import the stars feedback loop functions. */
#define FUNCTION feedback
Loic Hausammann's avatar
Loic Hausammann committed
125
#define FUNCTION_TASK_LOOP TASK_LOOP_FEEDBACK
Loic Hausammann's avatar
Loic Hausammann committed
126
#include "runner_doiact_stars.h"
Loic Hausammann's avatar
Loic Hausammann committed
127
#undef FUNCTION_TASK_LOOP
128
#undef FUNCTION
129

130
131
132
133
134
135
136
137
138
139
140
141
142
143
/* Import the black hole density loop functions. */
#define FUNCTION density
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
#include "runner_doiact_black_holes.h"
#undef FUNCTION_TASK_LOOP
#undef FUNCTION

/* Import the black hole feedback loop functions. */
#define FUNCTION feedback
#define FUNCTION_TASK_LOOP TASK_LOOP_FEEDBACK
#include "runner_doiact_black_holes.h"
#undef FUNCTION_TASK_LOOP
#undef FUNCTION

144
145
146
147
148
149
150
151
/**
 * @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
152
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
153

154
  struct spart *restrict sparts = c->stars.parts;
155
  const struct engine *e = r->e;
156
157
  const struct unit_system *us = e->internal_units;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
158
  const struct cosmology *cosmo = e->cosmology;
159
  const struct feedback_props *feedback_props = e->feedback_props;
160
161
162
163
164
165
  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;
166
  int redo = 0, scount = 0;
167

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

171
172
173
  TIMER_TIC;

  /* Anything to do here? */
174
  if (c->stars.count == 0) return;
Loic Hausammann's avatar
Loic Hausammann committed
175
  if (!cell_is_active_stars(c, e)) return;
176
177
178

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

Matthieu Schaller's avatar
Matthieu Schaller committed
183
184
        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->stars.h_max);
185
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
186
    }
187
188
189
190
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
191
192
193
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
194
    if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
Loic Hausammann's avatar
Loic Hausammann committed
195
      error("Can't allocate memory for sid.");
196
    if ((h_0 = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
197
      error("Can't allocate memory for h_0.");
198
    if ((left = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
199
      error("Can't allocate memory for left.");
200
    if ((right = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
201
      error("Can't allocate memory for right.");
202
    for (int k = 0; k < c->stars.count; k++)
203
204
      if (spart_is_active(&sparts[k], e) &&
          feedback_is_active(&sparts[k], e->time, cosmo, with_cosmology)) {
205
        sid[scount] = k;
206
207
208
        h_0[scount] = sparts[k].h;
        left[scount] = 0.f;
        right[scount] = stars_h_max;
209
        ++scount;
210
211
212
      }

    /* While there are particles that need to be updated... */
213
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
214
215
216
217
218
219
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
220
      for (int i = 0; i < scount; i++) {
221
222
223
224
225
226

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
227
228
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
229
230
231
#endif

        /* Get some useful values */
232
        const float h_init = h_0[i];
233
234
235
        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);
236

237
238
239
        float h_new;
        int has_no_neighbours = 0;

240
        if (sp->density.wcount == 0.f) { /* No neighbours case */
241
242
243
244
245
246

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

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

248
249
250
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
251
          stars_end_density(sp, cosmo);
252
253

          /* Compute one step of the Newton-Raphson scheme */
254
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
255
          const float n_target = stars_eta_dim;
256
257
          const float f = n_sum - n_target;
          const float f_prime =
258
259
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
260

261
          /* Improve the bisection bounds */
262
263
264
265
          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);
266
267
268
269
270
271
272

#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

273
          /* Skip if h is already h_max and we don't have enough neighbours */
274
275
276
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
277

278
            stars_reset_feedback(sp);
279
            feedback_reset_feedback(sp, feedback_props);
280

281
282
283
284
285
            /* Ok, we are done with this particle */
            continue;
          }

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

287
288
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
289
290
291
292
293
294
295
296
297
298
299
300

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

301
302
303
          /* 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);
304
305
306
307

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
308
309
310
311
312
313
        }

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

          /* Ok, correct then */
314
315
316
317
318
319
320
321
322
323
324
325
326
327

          /* 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;
          }
328
329

          /* If below the absolute maximum, try again */
330
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
331
332
333

            /* Flag for another round of fun */
            sid[redo] = sid[i];
334
335
336
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
337
338
339
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
340
            stars_init_spart(sp);
341
            feedback_init_spart(sp);
342
343
344
345

            /* Off we go ! */
            continue;

346
347
348
349
350
351
          } 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) {
352
353

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
354
            sp->h = stars_h_max;
355
356
357

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
358
              stars_spart_has_no_neighbours(sp, cosmo);
359
            }
360
361
362
363
364

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
365
366
367
          }
        }

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

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

373
        stars_reset_feedback(sp);
374

375
        /* Only do feedback if stars have a reasonable birth time */
376
        if (feedback_do_feedback(sp)) {
377

378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
          const integertime_t ti_step = get_integer_timestep(sp->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(e->ti_current - 1, sp->time_bin);

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

          /* Calculate age of the star at current time */
          double star_age_end_of_step;
393
          if (with_cosmology) {
394
            star_age_end_of_step = cosmology_get_delta_time_from_scale_factors(
395
396
                cosmo, sp->birth_scale_factor, (float)cosmo->a);
          } else {
397
            star_age_end_of_step = (float)e->time - sp->birth_time;
398
399
          }

400
401
402
403
404
405
406
407
408
409
          /* Has this star been around for a while ? */
          if (star_age_end_of_step > 0.) {

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

            /* Compute the stellar evolution  */
            feedback_evolve_spart(sp, feedback_props, cosmo, us,
                                  star_age_beg_of_step, dt);
410
411
412
413
          } else {

            /* Reset the feedback fields of the star particle */
            feedback_reset_feedback(sp, feedback_props);
414
          }
415
416
417
418
        } else {

          /* Reset the feedback fields of the star particle */
          feedback_reset_feedback(sp, feedback_props);
419
        }
420
421
422
423
424
425
      }

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

      /* Re-set the counter for the next loop (potentially). */
426
427
      scount = redo;
      if (scount > 0) {
428
429
430
431
432

        /* 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
433
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
434
435
436
437
438
439
440
441

#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)
442
443
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
444
445
446
447
448
449

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

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

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
459
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
460
                                                NULL, 1);
461
462
463
464
465
466

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

              /* Left or right? */
              if (l->t->ci == finger)
467
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
468
                                                  scount, l->t->cj, 1);
469
              else
470
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
471
                                                  scount, l->t->ci, 1);
472
473
474
475
476
477
            }
          }
        }
      }
    }

478
479
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
480
481
482
    }

    /* Be clean */
483
484
    free(left);
    free(right);
485
    free(sid);
486
    free(h_0);
487
488
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
489
490
  /* Update h_max */
  c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
491

492
  /* The ghost may not always be at the top level.
493
   * Therefore we need to update h_max between the super- and top-levels */
494
  if (c->stars.ghost) {
495
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
496
      atomic_max_d(&tmp->stars.h_max, h_max);
497
498
499
    }
  }

500
  if (timer) TIMER_TOC(timer_do_stars_ghost);
501
502
}

Tom Theuns's avatar
Tom Theuns committed
503
504
505
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
506
507
508
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
509
 */
510
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
511

512
513
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
514
515
516
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
517
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
518

519
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
520

521
  /* Anything to do here? */
522
  if (!cell_is_active_gravity(c, e)) return;
523

Tom Theuns's avatar
Tom Theuns committed
524
525
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
526
    for (int k = 0; k < 8; k++)
527
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
528
  } else {
529

530
531
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
532

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

536
      /* Is this part within the time step? */
537
      if (gpart_is_active(gp, e)) {
538
539
        external_gravity_acceleration(time, potential, constants, gp);
      }
540
    }
541
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
542

543
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
544
545
}

546
547
548
549
550
551
552
553
554
/**
 * @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) {

555
556
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
  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
581
/**
582
583
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
584
585
586
587
588
589
590
 *
 * @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) {

591
  const struct engine *e = r->e;
592
593
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
594
595
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
596
  const struct unit_system *us = e->internal_units;
597
  const struct hydro_props *hydro_props = e->hydro_properties;
598
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
599
  const double time_base = e->time_base;
600
  const integertime_t ti_current = e->ti_current;
601
602
603
  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
604
605
606

  TIMER_TIC;

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

Stefan Arridge's avatar
Stefan Arridge committed
610
611
612
613
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
614
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
615

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

619
620
621
      /* 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
622

623
      if (part_is_active(p, e)) {
624

625
        double dt_cool, dt_therm;
626
627
628
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
629
630
              get_integer_time_begin(ti_current - 1, p->time_bin);

631
632
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
633
634
635
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

636
637
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
638
          dt_therm = get_timestep(p->time_bin, time_base);
639
        }
640

641
        /* Let's cool ! */
642
643
644
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
645
      }
Stefan Arridge's avatar
Stefan Arridge committed
646
647
648
649
650
651
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
652
653
654
655
656
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

657
  struct engine *e = r->e;
658
  const struct cosmology *cosmo = e->cosmology;
659
660
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
661
662
663
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
664
  const int with_cosmology = (e->policy & engine_policy_cosmology);
665
  const int with_feedback = (e->policy & engine_policy_feedback);
666
667
668
  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;
669
  const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
670
671
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
672
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
673
674
675

  TIMER_TIC;

676
677
678
679
680
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID != e->nodeID)
    error("Running star formation task on a foreign node!");
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
681
  /* Anything to do here? */
682
  if (c->hydro.count == 0 || !cell_is_active_hydro(c, e)) {
683
    star_formation_logger_log_inactive_cell(&c->stars.sfh);
684
685
    return;
  }
686
687

  /* Reset the SFR */
688
  star_formation_logger_init(&c->stars.sfh);
Matthieu Schaller's avatar
Matthieu Schaller committed
689
690
691
692

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
693
694
695
696
697
698
699
700
      if (c->progeny[k] != NULL) {
        /* Load the child cell */
        struct cell *restrict cp = c->progeny[k];

        /* Do the recursion */
        runner_do_star_formation(r, cp, 0);

        /* Update current cell using child cells */
701
        star_formation_logger_add(&c->stars.sfh, &cp->stars.sfh);
702
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
703
  } else {
704
705
706
707
708
709
710
711

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

712
      /* Only work on active particles */
713
714
      if (part_is_active(p, e)) {

715
716
        /* 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
717
718
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
719

720
721
722
723
724
725
          /* 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);
726

727
728
729
730
731
732
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

734
735
736
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
737

738
          /* Add the SFR and SFR*dt to the SFH struct of this cell */
Folkert Nobels's avatar
Folkert Nobels committed
739
          star_formation_logger_log_active_part(p, xp, &c->stars.sfh, dt_star);
740

741
742
743
744
745
746
747
          /* 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);

748
749
750
751
752
            /* Did we get a star? (Or did we run out of spare ones?) */
            if (sp != NULL) {

              /* Copy the properties of the gas particle to the star particle */
              star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
Folkert Nobels's avatar
Folkert Nobels committed
753
754
                                             with_cosmology, phys_const,
                                             hydro_props, us, cooling);
755
756
757

              /* Update the Star formation history */
              star_formation_logger_log_new_spart(sp, &c->stars.sfh);
758
            }
759
760
761
762
763
764
765
766
767
          }

        } 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? */
768

Folkert Nobels's avatar
Folkert Nobels committed
769
      } else { /* is active? */
770

771
        /* Check if the particle is not inhibited */
772
        if (!part_is_inhibited(p, e)) {
773
774
          star_formation_logger_log_inactive_part(p, xp, &c->stars.sfh);
        }
775
      }
Folkert Nobels's avatar
Folkert Nobels committed
776
    } /* Loop over particles */
Matthieu Schaller's avatar
Matthieu Schaller committed
777
778
  }

779
780
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
781
  if (with_feedback && (c == c->top) &&
782
      (current_stars_count != c->stars.count)) {
783
784

    cell_clear_stars_sort_flags(c);
785
    runner_do_all_stars_sort(r, c);
786
  }
787

Matthieu Schaller's avatar
Matthieu Schaller committed
788
789
790
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
791
792
793
794
795
796
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
797
void runner_do_sort_ascending(struct entry *sort, int N) {
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851

  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
852
        }
853
854
855
856
857
858
859
860
861
862
863
864
      } 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
865
    }
866
867
868
  }
}

869
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
870
871
872
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
873
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
874
 */
875
876
877
878
879
880
881
882
883
#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
884
#else
Loic Hausammann's avatar
Loic Hausammann committed
885
886
887
888
#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
889
#endif
Loic Hausammann's avatar
Loic Hausammann committed
890
891
892

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
893

Pedro Gonnet's avatar
Pedro Gonnet committed
894
895
896
897
898
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
899
 * @param flags Cell flag.
900
901
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
902
903
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
904
 */
Loic Hausammann's avatar
Loic Hausammann committed
905
906
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
907
908

  struct entry *fingers[8];
909
910
911
  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
912
  float buff[8];
913

914
915
  TIMER_TIC;

916
917
918
919
#ifdef SWIFT_DEBUG_CHECKS
  if (c->hydro.super == NULL) error("Task called above the super level!!!");
#endif

920
  /* We need to do the local sorts plus whatever was requested further up. */
921
  flags |= c->hydro.do_sort;
922
  if (cleanup) {
923
    c->hydro.sorted = 0;
924
  } else {
925
    flags &= ~c->hydro.sorted;
926
  }
927
  if (flags == 0 && !cell_get_flag(c, cell_flag_do_hydro_sub_sort)) return;
928
929

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

933
934
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
935
  runner_check_sorts_hydro(c, c->hydro.sorted);
936
937

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
938
939
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
940
941
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
942
  }
943

944
945
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
946
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
947
#endif
948

949
950
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
951
952
953
954
955

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

    /* Fill in the gaps within the progeny. */
956
    float dx_max_sort = 0.0f;
957
    float dx_max_sort_old = 0.0f;
958
    for (int k = 0; k < 8; k++) {
959
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
960
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
961
        runner_do_hydro_sort(r, c->progeny[k], flags,
962
963
964
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
965
966
967
        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);
968
      }
Pedro Gonnet's avatar