runner.c 119 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 "star_formation_logger.h"
68
#include "stars.h"
69
70
#include "task.h"
#include "timers.h"
71
#include "timestep.h"
72
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
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
140
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
141
142
143
144
145
146
  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;
147
  int redo = 0, scount = 0;
148

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

152
153
154
  TIMER_TIC;

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

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

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

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

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

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

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

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
207
208
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
209
210
211
#endif

        /* Get some useful values */
212
        const float h_init = h_0[i];
213
214
215
        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);
216

217
218
219
        float h_new;
        int has_no_neighbours = 0;

220
        if (sp->density.wcount == 0.f) { /* No neighbours case */
221
222
223
224
225
226

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

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

228
229
230
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
231
          stars_end_density(sp, cosmo);
232
233

          /* Compute one step of the Newton-Raphson scheme */
234
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
235
          const float n_target = stars_eta_dim;
236
237
          const float f = n_sum - n_target;
          const float f_prime =
238
239
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
240

241
          /* Improve the bisection bounds */
242
243
244
245
          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);
246
247
248
249
250
251
252

#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

253
          /* Skip if h is already h_max and we don't have enough neighbours */
254
255
256
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
257

258
259
            stars_reset_feedback(sp);

260
261
262
263
264
            /* Ok, we are done with this particle */
            continue;
          }

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

266
267
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
268
269
270
271
272
273
274
275
276
277
278
279

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

280
281
282
283
284
285
286
287
288
#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);
289
290
291
292

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
293
294
295
296
297
298
        }

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

          /* Ok, correct then */
299
300
301
302
303
304
305
306
307
308
309
310
311
312

          /* 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;
          }
313
314

          /* If below the absolute maximum, try again */
315
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
316
317
318

            /* Flag for another round of fun */
            sid[redo] = sid[i];
319
320
321
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
322
323
324
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
325
            stars_init_spart(sp);
326
327
328
329

            /* Off we go ! */
            continue;

330
331
332
333
334
335
          } 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) {
336
337

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
338
            sp->h = stars_h_max;
339
340
341

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
342
              stars_spart_has_no_neighbours(sp, cosmo);
343
            }
344
345
346
347
348

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
349
350
351
          }
        }

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

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

357
        stars_reset_feedback(sp);
358

Loic Hausammann's avatar
Loic Hausammann committed
359
        /* Compute the stellar evolution  */
360
        stars_evolve_spart(sp, e->stars_properties, cosmo);
361
362
363
364
365
366
      }

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

      /* Re-set the counter for the next loop (potentially). */
367
368
      scount = redo;
      if (scount > 0) {
369
370
371
372
373

        /* 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
374
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
375
376
377
378
379
380
381
382

#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)
383
384
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
385
386
387
388
389
390

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

              /* Left or right? */
              if (l->t->ci == finger)
391
392
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
393
              else
394
395
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
396
397
398
399
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
400
401
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
402
403
404
405
406
407

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

              /* Left or right? */
              if (l->t->ci == finger)
408
409
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
410
              else
411
412
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
413
414
415
416
417
418
            }
          }
        }
      }
    }

419
420
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
421
422
423
    }

    /* Be clean */
424
425
    free(left);
    free(right);
426
    free(sid);
427
    free(h_0);
428
429
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
430
431
  /* Update h_max */
  c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
432

433
  /* The ghost may not always be at the top level.
434
   * Therefore we need to update h_max between the super- and top-levels */
435
  if (c->stars.ghost) {
436
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
437
      atomic_max_d(&tmp->stars.h_max, h_max);
438
439
440
    }
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
441
  if (timer) TIMER_TOC(timer_dostars_ghost);
442
443
}

Tom Theuns's avatar
Tom Theuns committed
444
445
446
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
447
448
449
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
450
 */
451
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
452

453
454
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
455
456
457
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
458
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
459

460
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
461

462
  /* Anything to do here? */
463
  if (!cell_is_active_gravity(c, e)) return;
464

Tom Theuns's avatar
Tom Theuns committed
465
466
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
467
    for (int k = 0; k < 8; k++)
468
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
469
  } else {
470

471
472
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
473

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

477
      /* Is this part within the time step? */
478
      if (gpart_is_active(gp, e)) {
479
480
        external_gravity_acceleration(time, potential, constants, gp);
      }
481
    }
482
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
483

484
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
485
486
}

487
488
489
490
491
492
493
494
495
/**
 * @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) {

496
497
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
  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
522
/**
523
524
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
525
526
527
528
529
530
531
 *
 * @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) {

532
  const struct engine *e = r->e;
533
534
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
535
536
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
537
  const struct unit_system *us = e->internal_units;
538
  const struct hydro_props *hydro_props = e->hydro_properties;
539
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
540
  const double time_base = e->time_base;
541
  const integertime_t ti_current = e->ti_current;
542
543
544
  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
545
546
547

  TIMER_TIC;

548
  /* Anything to do here? */
549
  if (!cell_is_active_hydro(c, e)) return;
550

Stefan Arridge's avatar
Stefan Arridge committed
551
552
553
554
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
555
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
556

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

560
561
562
      /* 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
563

564
      if (part_is_active(p, e)) {
565

566
        double dt_cool, dt_therm;
567
568
569
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
570
571
              get_integer_time_begin(ti_current - 1, p->time_bin);

572
573
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
574
575
576
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

577
578
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
579
          dt_therm = get_timestep(p->time_bin, time_base);
580
        }
581

582
        /* Let's cool ! */
583
584
585
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
586
      }
Stefan Arridge's avatar
Stefan Arridge committed
587
588
589
590
591
592
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
593
594
595
596
597
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

598
  struct engine *e = r->e;
599
  const struct cosmology *cosmo = e->cosmology;
600
601
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
602
603
604
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
605
  const int with_cosmology = (e->policy & engine_policy_cosmology);
606
  const int with_feedback = (e->policy & engine_policy_feedback);
607
608
609
  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;
610
  const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
611
612
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
613
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
614
615
616

  TIMER_TIC;

617
618
619
620
621
#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
622
  /* Anything to do here? */
623
  if (c->hydro.count == 0) return;
624
  if (!cell_is_active_hydro(c, e)) {
625
    star_formation_logger_log_inactive_cell(&c->stars.sfh);
626
627
    return;
  }
628
  star_formation_logger_init(&c->stars.sfh);
Matthieu Schaller's avatar
Matthieu Schaller committed
629
630
631
632

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
633
634
635
636
637
638
639
640
      if (c->progeny[k] != NULL) {
        /* Load the child cell */
        struct cell *restrict cp = c->progeny[k];

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

        /* Update current cell using child cells */
Folkert Nobels's avatar
Folkert Nobels committed
641
642
        star_formation_logger_add_first_to_second(&cp->stars.sfh,
                                                  &c->stars.sfh);
643
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
644
  } else {
645
646
647
648
649
650
651
652

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

653
      /* Only work on active particles */
654
655
      if (part_is_active(p, e)) {

656
657
        /* 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
658
659
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
660

661
662
663
664
665
666
          /* 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);
667

668
669
670
671
672
673
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

675
676
677
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
678

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

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

689
690
691
692
693
694
            /* 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);
695
696
697

              /* Update the Star formation history */
              star_formation_logger_log_new_spart(sp, &c->stars.sfh);
698
            }
699
700
701
702
703
704
705
706
          }

        } 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);

Folkert Nobels's avatar
Folkert Nobels committed
707
708
        }      /* Not Star-forming? */
      } else { /* is active? */
709
        /* Check if the particle is not inhibited */
710
        if (!part_is_inhibited(p, e)) {
711
712
          star_formation_logger_log_inactive_part(p, xp, &c->stars.sfh);
        }
713
      }
Folkert Nobels's avatar
Folkert Nobels committed
714
    } /* Loop over particles */
Matthieu Schaller's avatar
Matthieu Schaller committed
715
716
  }

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

    cell_clear_stars_sort_flags(c);
723
    runner_do_all_stars_sort(r, c);
724
  }
725

Matthieu Schaller's avatar
Matthieu Schaller committed
726
727
728
  if (timer) TIMER_TOC(timer_do_star_formation);
}

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

  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
790
        }
791
792
793
794
795
796
797
798
799
800
801
802
      } 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
803
    }
804
805
806
  }
}

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

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
831

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

  struct entry *fingers[8];
847
848
849
  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
850
  float buff[8];
851

852
853
  TIMER_TIC;

854
855
856
857
#ifdef SWIFT_DEBUG_CHECKS
  if (c->hydro.super == NULL) error("Task called above the super level!!!");
#endif

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

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

871
872
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
873
  runner_check_sorts_hydro(c, c->hydro.sorted);
874
875

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

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

887
888
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
889
890
891
892
893

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

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

    /* Loop over the 13 different sort arrays. */
912
    for (int j = 0; j < 13; j++) {
913
914
915
916
917

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

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

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

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

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

        /* 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. */
960
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
961
          int temp_i = inds[k - 1];
962
963
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar