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
#include "feedback.h"
54
55
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
56
#include "hydro_properties.h"
57
#include "kick.h"
Loikki's avatar
Loikki committed
58
#include "logger.h"
59
#include "memuse.h"
60
#include "minmax.h"
James Willis's avatar
James Willis committed
61
#include "runner_doiact_vec.h"
62
#include "scheduler.h"
63
#include "sort_part.h"
64
#include "space.h"
65
#include "space_getsid.h"
66
#include "star_formation.h"
67
#include "star_formation_iact.h"
68
#include "stars.h"
69
70
#include "task.h"
#include "timers.h"
71
#include "timestep.h"
72
#include "timestep_limiter.h"
73
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
74

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

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

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

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

104
105
106
107
108
109
110
/* 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

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

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

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

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

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

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

155
156
157
  TIMER_TIC;

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

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

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

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

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

199
200
201
202
      /* Reset the redo-count. */
      redo = 0;

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

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

208
209
210
211
212
213
214
215
216
217
218
219
        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);
        }
Alexei Borissov's avatar
Alexei Borissov committed
220

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

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

233
234
235
        float h_new;
        int has_no_neighbours = 0;

236
        if (sp->density.wcount == 0.f) { /* No neighbours case */
237
238
239
240
241
242

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

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

244
245
246
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
247
          stars_end_density(sp, cosmo);
248
249

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

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

#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

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

274
            stars_reset_feedback(sp);
275
276
277
278
279
280

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

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

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

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

296
297
298
299
300
301
302
303
304
#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);
305
306
307
308

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

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

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

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

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

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

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
341
            stars_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
377
378
379
380
381
382
383
384
385
386
387
388
389
        if (sp->birth_time != -1.) {

          /* Calculate age of the star */
          double star_age;
          if (with_cosmology) {
            star_age = cosmology_get_delta_time_from_scale_factors(
                cosmo, sp->birth_scale_factor, (float)cosmo->a);
          } else {
            star_age = e->time - sp->birth_time;
          }

          /* Compute the stellar evolution  */
          feedback_evolve_spart(sp, feedback_props, cosmo, us, star_age, dt);
        }
390
391
392
393
394
395
      }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

701
702
703
704
705
706
707
            /* 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);
            }
708
709
710
711
712
713
714
715
716
717
718
          }

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

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

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

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

Pedro Gonnet's avatar
Pedro Gonnet committed
733
734
735
736
737
738
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
739
void runner_do_sort_ascending(struct entry *sort, int N) {
740
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

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

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

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
835

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

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

856
857
  TIMER_TIC;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        /* 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. */
964
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
965
          int temp_i = inds[k - 1];
966
967
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
968
        }
969

970
971
972
      } /* Merge. */

      /* Add a sentinel. */
973
974
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
975
976

      /* Mark as sorted. */
977