runner.c 116 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
140
  // const struct unit_system *us = e->internal_units;
  // 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 */
207
208
209
210
211
212
213
214
215
216
217
        /* double dt; */
        /* if (with_cosmology) { */
        /*   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); */
        /*   dt = cosmology_get_therm_kick_factor(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
218

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

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

231
232
233
        float h_new;
        int has_no_neighbours = 0;

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

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

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

242
243
244
        } else {

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

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

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

#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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

            /* Off we go ! */
            continue;

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

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

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

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

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

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

371
        stars_reset_feedback(sp);
372

373
        // MATTHIEU
374

375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
        /* Only do feedback if stars have a reasonable birth time */
        /* if (sp->birth_time != -1) { */

        /*   /\* Calculate age of the star *\/ */
        /*   double star_age, current_time_begin = -1; */
        /*   if (with_cosmology) { */
        /*     star_age = cosmology_get_delta_time_from_scale_factors( */
        /*         cosmo, sp->birth_scale_factor, (float)cosmo->a); */
        /*   } else { */
        /*     current_time_begin = */
        /*         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; */
        /*   } */

        /*   /\* Compute the stellar evolution  *\/ */
        /*   //stars_evolve_spart(sp, e->stars_properties, cosmo, us, star_age,
         * dt); */
        /* } */
395
396
397
398
399
400
      }

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

      /* Re-set the counter for the next loop (potentially). */
401
402
      scount = redo;
      if (scount > 0) {
403
404
405
406
407

        /* 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
408
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
409
410
411
412
413
414
415
416

#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)
417
418
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
419
420
421
422
423
424

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

              /* Left or right? */
              if (l->t->ci == finger)
425
426
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
427
              else
428
429
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
430
431
432
433
            }

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

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

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

453
454
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
455
456
457
    }

    /* Be clean */
458
459
    free(left);
    free(right);
460
    free(sid);
461
    free(h_0);
462
463
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
464
465
  /* Update h_max */
  c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
466

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

Matthieu Schaller's avatar
Matthieu Schaller committed
475
  if (timer) TIMER_TOC(timer_dostars_ghost);
476
477
}

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

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

494
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
495

496
  /* Anything to do here? */
497
  if (!cell_is_active_gravity(c, e)) return;
498

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

505
506
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
507

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

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

518
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
519
520
}

521
522
523
524
525
526
527
528
529
/**
 * @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) {

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

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

  TIMER_TIC;

582
  /* Anything to do here? */
583
  if (!cell_is_active_hydro(c, e)) return;
584

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

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

594
595
596
      /* 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
597

598
      if (part_is_active(p, e)) {
599

600
        double dt_cool, dt_therm;
601
602
603
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
604
605
              get_integer_time_begin(ti_current - 1, p->time_bin);

606
607
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
608
609
610
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

611
612
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
613
          dt_therm = get_timestep(p->time_bin, time_base);
614
        }
615

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
627
628
629
630
631
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

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

  TIMER_TIC;

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

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

673
      /* Only work on active particles */
674
675
      if (part_is_active(p, e)) {

676
677
        /* 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
678
679
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
680

681
682
683
684
685
686
          /* 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);
687

688
689
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
690

691
692
693
          } else {
            dt_star = get_timestep(p->time_bin, time_base);
          }
694

695
696
697
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
698

699
700
701
702
703
704
705
          /* 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);

706
707
708
709
710
711
712
            /* 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);
            }
713
714
715
716
717
718
719
720
721
722
723
          }

        } 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
724
725
  }

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

    cell_clear_stars_sort_flags(c);
732
    runner_do_all_stars_sort(r, c);
733
  }
734

Matthieu Schaller's avatar
Matthieu Schaller committed
735
736
737
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
738
739
740
741
742
743
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
744
void runner_do_sort_ascending(struct entry *sort, int N) {
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
795
796
797
798

  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
799
        }
800
801
802
803
804
805
806
807
808
809
810
811
      } 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
812
    }
813
814
815
  }
}

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

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
840

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

  struct entry *fingers[8];
856
857
858
  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
859
  float buff[8];
860

861
862
  TIMER_TIC;

863
864
865
866
#ifdef SWIFT_DEBUG_CHECKS
  if (c->hydro.super == NULL) error("Task called above the super level!!!");
#endif

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

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

880
881
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
882
  runner_check_sorts_hydro(c, c->hydro.sorted);
883
884

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

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

896
897
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
898
899
900
901
902

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

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

    /* Loop over the 13 different sort arrays. */
921
    for (int j = 0; j < 13; j++) {
922
923
924
925
926

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

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

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

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

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

        /* 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. */
969
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
970
          int temp_i = inds[k - 1];
971
972
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
973
        }
974

975
976
977
      } /* Merge. */

      /* Add a sentinel. */
978
979
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
980
981

      /* Mark as sorted. */
982
      atomic_or(&c->hydro.sorted, 1 << j);