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
#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
144
145
146
147
148
  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;
149
  int redo = 0, scount = 0;
150

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

154
155
156
  TIMER_TIC;

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

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

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

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

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

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

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

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

207
        /* get particle timestep */
208
209
210
211
212
213
214
215
216
217
218
        /* 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
219

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

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

232
233
234
        float h_new;
        int has_no_neighbours = 0;

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

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

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

243
244
245
        } else {

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

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

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

#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

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

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

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

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

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

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

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

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

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

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

          /* Case where we have been oscillating around the solution */
          if ((h_new == left[i] && h_old == right[i]) ||
              (h_old == left[i] && h_new == right[i])) {

            /* Bissect the remaining interval */
            sp->h = pow_inv_dimension(
                0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));

          } else {

            /* Normal case */
            sp->h = h_new;
          }
328
329

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

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

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

            /* Off we go ! */
            continue;

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

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

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

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

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

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

372
        stars_reset_feedback(sp);
373

374
        // MATTHIEU
375

376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
        /* 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); */
        /* } */
396
397
398
399
400
401
      }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
841

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

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

862
863
  TIMER_TIC;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

976
977
978
      } /* Merge. */

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

      /* Mark as sorted. */