runner.c 115 KB
Newer Older
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
5
6
7
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
 *               2016 John A. Regan (john.a.regan@durham.ac.uk)
 *                    Tom Theuns (tom.theuns@durham.ac.uk)
8
 *
9
10
11
12
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
13
 *
14
15
16
17
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
18
 *
19
20
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
21
 *
22
 ******************************************************************************/
Pedro Gonnet's avatar
Pedro Gonnet committed
23

Pedro Gonnet's avatar
Pedro Gonnet committed
24
25
/* Config parameters. */
#include "../config.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
26
27
28
29

/* Some standard headers. */
#include <float.h>
#include <limits.h>
30
#include <stdlib.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
31

32
33
/* MPI headers. */
#ifdef WITH_MPI
34
#include <mpi.h>
35
36
#endif

37
38
39
/* This object's header. */
#include "runner.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
40
/* Local headers. */
41
#include "active.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
42
#include "approx_math.h"
43
#include "atomic.h"
44
#include "cell.h"
45
#include "chemistry.h"
46
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
47
#include "cooling.h"
48
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
49
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
50
#include "engine.h"
51
#include "entropy_floor.h"
52
#include "error.h"
53
54
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
55
#include "hydro_properties.h"
56
#include "kick.h"
Loikki's avatar
Loikki committed
57
#include "logger.h"
58
#include "memuse.h"
59
#include "minmax.h"
James Willis's avatar
James Willis committed
60
#include "runner_doiact_vec.h"
61
#include "scheduler.h"
62
#include "sort_part.h"
63
#include "space.h"
64
#include "space_getsid.h"
65
#include "star_formation.h"
66
#include "star_formation_iact.h"
67
#include "stars.h"
68
69
#include "task.h"
#include "timers.h"
70
#include "timestep.h"
71
#include "timestep_limiter.h"
72
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
73

74
75
76
77
#define TASK_LOOP_DENSITY 0
#define TASK_LOOP_GRADIENT 1
#define TASK_LOOP_FORCE 2
#define TASK_LOOP_LIMITER 3
Loic Hausammann's avatar
Loic Hausammann committed
78
#define TASK_LOOP_FEEDBACK 4
79

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

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

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

103
104
105
106
107
108
109
/* Import the limiter loop functions. */
#define FUNCTION limiter
#define FUNCTION_TASK_LOOP TASK_LOOP_LIMITER
#include "runner_doiact.h"
#undef FUNCTION
#undef FUNCTION_TASK_LOOP

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

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

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

/**
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
Loic Hausammann's avatar
Loic Hausammann committed
135
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
136

137
  struct spart *restrict sparts = c->stars.parts;
138
  const struct engine *e = r->e;
139
  const struct unit_system *us = e->internal_units;
Alexei Borissov's avatar
Alexei Borissov committed
140
  const int with_cosmology = (e->policy & engine_policy_cosmology);
141
  const struct cosmology *cosmo = e->cosmology;
142
143
144
145
146
147
  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;
148
  int redo = 0, scount = 0;
149

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

153
154
155
  TIMER_TIC;

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
165
166
        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->stars.h_max);
167
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
168
    }
169
170
171
172
  } else {

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

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

197
198
199
200
      /* Reset the redo-count. */
      redo = 0;

      /* Loop over the remaining active parts in this cell. */
201
      for (int i = 0; i < scount; i++) {
202
203
204
205

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

206
        /* get particle timestep */
Alexei Borissov's avatar
Alexei Borissov committed
207
208
209
210
        double dt;
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(sp->time_bin);
          const integertime_t ti_begin =
Alexei Borissov's avatar
Alexei Borissov committed
211
              get_integer_time_begin(e->ti_current - 1, sp->time_bin);
Alexei Borissov's avatar
Alexei Borissov committed
212
          dt = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
213
                                               ti_begin + ti_step);
Alexei Borissov's avatar
Alexei Borissov committed
214
        } else {
Alexei Borissov's avatar
Alexei Borissov committed
215
          dt = get_timestep(sp->time_bin, e->time_base);
Alexei Borissov's avatar
Alexei Borissov committed
216
217
        }

218
219
#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
220
221
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
222
223
224
#endif

        /* Get some useful values */
225
        const float h_init = h_0[i];
226
227
228
        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);
229

230
231
232
        float h_new;
        int has_no_neighbours = 0;

233
        if (sp->density.wcount == 0.f) { /* No neighbours case */
234
235
236
237
238
239

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

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

241
242
243
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
244
          stars_end_density(sp, cosmo);
245
246

          /* Compute one step of the Newton-Raphson scheme */
247
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
248
          const float n_target = stars_eta_dim;
249
250
          const float f = n_sum - n_target;
          const float f_prime =
251
252
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
253

254
          /* Improve the bisection bounds */
255
256
257
258
          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);
259
260
261
262
263
264
265

#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

266
          /* Skip if h is already h_max and we don't have enough neighbours */
267
268
269
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
270

271
            stars_reset_feedback(sp);
272
273
274
275
276
277

            /* Ok, we are done with this particle */
            continue;
          }

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

279
280
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
281
282
283
284
285
286
287
288
289
290
291
292

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

293
294
295
296
297
298
299
300
301
#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);
302
303
304
305

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

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

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

          /* 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;
          }
326
327

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

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

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
338
            stars_init_spart(sp);
339
340
341
342

            /* Off we go ! */
            continue;

343
344
345
346
347
348
          } 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) {
349
350

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

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

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
362
363
364
          }
        }

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

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

370
        stars_reset_feedback(sp);
371

372
373
        /* Only do feedback if stars have a reasonable birth time */
        if (sp->birth_time != -1) {
374
375

          /* Calculate age of the star */
376
377
378
          double star_age, current_time_begin = -1;
          if (with_cosmology) {
            star_age = cosmology_get_delta_time_from_scale_factors(
379
                cosmo, sp->birth_scale_factor, (float)cosmo->a);
380
          } else {
381
            current_time_begin =
382
383
384
385
386
387
                get_integer_time_begin(e->ti_current - 1, sp->time_bin) *
                    e->time_base +
                e->time_begin;
            star_age = current_time_begin - sp->birth_time;
          }

388
          /* Compute the stellar evolution  */
389
          stars_evolve_spart(sp, e->stars_properties, cosmo, us, star_age, dt);
390
        }
391
392
393
394
395
396
      }

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

      /* Re-set the counter for the next loop (potentially). */
397
398
      scount = redo;
      if (scount > 0) {
399
400
401
402
403

        /* 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
404
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
405
406
407
408
409
410
411
412

#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)
413
414
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
415
416
417
418
419
420

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

              /* Left or right? */
              if (l->t->ci == finger)
421
422
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
423
              else
424
425
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
426
427
428
429
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
430
431
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
432
433
434
435
436
437

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

              /* Left or right? */
              if (l->t->ci == finger)
438
439
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
440
              else
441
442
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
443
444
445
446
447
448
            }
          }
        }
      }
    }

449
450
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
451
452
453
    }

    /* Be clean */
454
455
    free(left);
    free(right);
456
    free(sid);
457
    free(h_0);
458
459
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
460
461
  /* Update h_max */
  c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
462

463
  /* The ghost may not always be at the top level.
464
   * Therefore we need to update h_max between the super- and top-levels */
465
  if (c->stars.ghost) {
466
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
467
      atomic_max_d(&tmp->stars.h_max, h_max);
468
469
470
    }
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
471
  if (timer) TIMER_TOC(timer_dostars_ghost);
472
473
}

Tom Theuns's avatar
Tom Theuns committed
474
475
476
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
477
478
479
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
480
 */
481
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
482

483
484
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
485
486
487
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
488
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
489

490
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
491

492
  /* Anything to do here? */
493
  if (!cell_is_active_gravity(c, e)) return;
494

Tom Theuns's avatar
Tom Theuns committed
495
496
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
497
    for (int k = 0; k < 8; k++)
498
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
499
  } else {
500

501
502
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
503

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

507
      /* Is this part within the time step? */
508
      if (gpart_is_active(gp, e)) {
509
510
        external_gravity_acceleration(time, potential, constants, gp);
      }
511
    }
512
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
513

514
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
515
516
}

517
518
519
520
521
522
523
524
525
/**
 * @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) {

526
527
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
  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
552
/**
553
554
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
555
556
557
558
559
560
561
 *
 * @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) {

562
  const struct engine *e = r->e;
563
564
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
565
566
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
567
  const struct unit_system *us = e->internal_units;
568
  const struct hydro_props *hydro_props = e->hydro_properties;
569
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
570
  const double time_base = e->time_base;
571
  const integertime_t ti_current = e->ti_current;
572
573
574
  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
575
576
577

  TIMER_TIC;

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

Stefan Arridge's avatar
Stefan Arridge committed
581
582
583
584
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
585
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
586

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

590
591
592
      /* 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
593

594
      if (part_is_active(p, e)) {
595

596
        double dt_cool, dt_therm;
597
598
599
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
600
601
              get_integer_time_begin(ti_current - 1, p->time_bin);

602
603
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
604
605
606
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

607
608
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
609
          dt_therm = get_timestep(p->time_bin, time_base);
610
        }
611

612
        /* Let's cool ! */
613
614
615
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
616
      }
Stefan Arridge's avatar
Stefan Arridge committed
617
618
619
620
621
622
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
623
624
625
626
627
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

628
  struct engine *e = r->e;
629
  const struct cosmology *cosmo = e->cosmology;
630
631
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
632
633
634
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
635
  const int with_cosmology = (e->policy & engine_policy_cosmology);
636
  const int with_feedback = (e->policy & engine_policy_feedback);
637
638
639
  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;
640
  const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
641
642
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
643
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
644
645
646

  TIMER_TIC;

647
648
649
650
651
#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
652
  /* Anything to do here? */
653
  if (c->hydro.count == 0) return;
Matthieu Schaller's avatar
Matthieu Schaller committed
654
655
656
657
658
659
660
  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 {
661
662
663
664
665
666
667
668

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

669
      /* Only work on active particles */
670
671
      if (part_is_active(p, e)) {

672
673
        /* 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
674
675
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
676

677
678
679
680
681
682
          /* 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);
683

684
685
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
686

687
688
689
          } else {
            dt_star = get_timestep(p->time_bin, time_base);
          }
690

691
692
693
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
694

695
696
697
698
699
700
701
          /* 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);

702
703
704
705
706
707
708
            /* Did we get a star? (Or did we run out of spare ones?) */
            if (sp != NULL) {

              /* Copy the properties of the gas particle to the star particle */
              star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
                                             with_cosmology);
            }
709
710
711
712
713
714
715
716
717
718
719
          }

        } 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
720
721
  }

722
723
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
724
  if (with_feedback && (c == c->top) &&
725
      (current_stars_count != c->stars.count)) {
726
727

    cell_clear_stars_sort_flags(c);
728
    runner_do_all_stars_sort(r, c);
729
  }
730

Matthieu Schaller's avatar
Matthieu Schaller committed
731
732
733
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
734
735
736
737
738
739
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
740
void runner_do_sort_ascending(struct entry *sort, int N) {
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794

  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
795
        }
796
797
798
799
800
801
802
803
804
805
806
807
      } 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
808
    }
809
810
811
  }
}

812
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
813
814
815
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
816
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
817
 */
818
819
820
821
822
823
824
825
826
#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
827
#else
Loic Hausammann's avatar
Loic Hausammann committed
828
829
830
831
#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
832
#endif
Loic Hausammann's avatar
Loic Hausammann committed
833
834
835

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
836

Pedro Gonnet's avatar
Pedro Gonnet committed
837
838
839
840
841
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
842
 * @param flags Cell flag.
843
844
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
845
846
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
847
 */
Loic Hausammann's avatar
Loic Hausammann committed
848
849
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
850
851

  struct entry *fingers[8];
852
853
854
  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
855
  float buff[8];
856

857
858
  TIMER_TIC;

859
860
861
862
#ifdef SWIFT_DEBUG_CHECKS
  if (c->hydro.super == NULL) error("Task called above the super level!!!");
#endif

863
  /* We need to do the local sorts plus whatever was requested further up. */
864
  flags |= c->hydro.do_sort;
865
  if (cleanup) {
866
    c->hydro.sorted = 0;
867
  } else {
868
    flags &= ~c->hydro.sorted;
869
  }
870
  if (flags == 0 && !c->hydro.do_sub_sort) return;
871
872

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

876
877
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
878
  runner_check_sorts_hydro(c, c->hydro.sorted);
879
880

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
881
882
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
883
884
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
885
  }
886

887
888
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
889
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
890
#endif
891

892
893
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
894
895
896
897
898

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

    /* Fill in the gaps within the progeny. */
899
    float dx_max_sort = 0.0f;
900
    float dx_max_sort_old = 0.0f;
901
    for (int k = 0; k < 8; k++) {
902
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
903
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
904
        runner_do_hydro_sort(r, c->progeny[k], flags,
905
906
907
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
908
909
910
        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);
911
      }
912
    }
913
914
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
915
916

    /* Loop over the 13 different sort arrays. */
917
    for (int j = 0; j < 13; j++) {
918
919
920
921
922

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

      /* Init the particle index offsets. */
923
      int off[8];
924
925
      off[0] = 0;
      for (int k = 1; k < 8; k++)
926
        if (c->progeny[k - 1] != NULL)
927
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
928
929
930
931
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
932
      int inds[8];
933
      for (int k = 0; k < 8; k++) {
934
        inds[k] = k;
935
936
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
937
938
939
940
941
942
943
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
944
945
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
946
          if (buff[inds[k]] < buff[inds[i]]) {
947
            int temp_i = inds[i];
948
949
950
951
952
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
953
      struct entry *finger = c->hydro.sort[j];
954
      for (int ind = 0; ind < count; ind++) {
955
956
957
958
959
960
961
962
963
964