runner.c 110 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
149
150
151

  TIMER_TIC;

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

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

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

    /* While there are particles that need to be updated... */
184
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
185
186
187
188
189
190
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
191
      for (int i = 0; i < scount; i++) {
192
193
194
195
196
197

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
198
199
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
200
201
202
#endif

        /* Get some useful values */
203
        const float h_init = h_0[i];
204
205
206
        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);
207

208
209
210
        float h_new;
        int has_no_neighbours = 0;

211
        if (sp->density.wcount == 0.f) { /* No neighbours case */
212
213
214
215
216
217

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

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

219
220
221
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
222
          stars_end_density(sp, cosmo);
223
224

          /* Compute one step of the Newton-Raphson scheme */
225
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
226
          const float n_target = stars_eta_dim;
227
228
          const float f = n_sum - n_target;
          const float f_prime =
229
230
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
231

232
          /* Improve the bisection bounds */
233
234
235
236
          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);
237
238
239
240
241
242
243

#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

244
          /* Skip if h is already h_max and we don't have enough neighbours */
245
246
247
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
248

249
250
            stars_reset_feedback(sp);

251
252
253
254
255
            /* Ok, we are done with this particle */
            continue;
          }

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

257
258
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
259
260
261
262
263
264
265
266
267
268
269
270

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

271
272
273
274
275
276
277
278
279
#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);
280
281
282
283

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
284
285
286
287
288
289
        }

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

          /* Ok, correct then */
290
291
292
293
294
295
296
297
298
299
300
301
302
303

          /* 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;
          }
304
305

          /* If below the absolute maximum, try again */
306
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
307
308
309

            /* Flag for another round of fun */
            sid[redo] = sid[i];
310
311
312
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
313
314
315
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
316
            stars_init_spart(sp);
317
318
319
320

            /* Off we go ! */
            continue;

321
322
323
324
325
326
          } 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) {
327
328

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
329
            sp->h = stars_h_max;
330
331
332

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
333
              stars_spart_has_no_neighbours(sp, cosmo);
334
            }
335
336
337
338
339

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
340
341
342
          }
        }

343
        /* We now have a particle whose smoothing length has converged */
344
        stars_reset_feedback(sp);
345

Loic Hausammann's avatar
Loic Hausammann committed
346
        /* Compute the stellar evolution  */
347
        stars_evolve_spart(sp, e->stars_properties, cosmo);
348
349
350
351
352
353
      }

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

      /* Re-set the counter for the next loop (potentially). */
354
355
      scount = redo;
      if (scount > 0) {
356
357
358
359
360

        /* 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
361
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
362
363
364
365
366
367
368
369

#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)
370
371
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
372
373
374
375
376
377

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

              /* Left or right? */
              if (l->t->ci == finger)
378
379
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
380
              else
381
382
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
383
384
385
386
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
387
388
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
389
390
391
392
393
394

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

              /* Left or right? */
              if (l->t->ci == finger)
395
396
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
397
              else
398
399
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
400
401
402
403
404
405
            }
          }
        }
      }
    }

406
407
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
408
409
410
    }

    /* Be clean */
411
412
    free(left);
    free(right);
413
    free(sid);
414
    free(h_0);
415
416
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
417
  if (timer) TIMER_TOC(timer_dostars_ghost);
418
419
}

Tom Theuns's avatar
Tom Theuns committed
420
421
422
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
423
424
425
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
426
 */
427
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
428

429
430
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
431
432
433
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
434
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
435

436
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
437

438
  /* Anything to do here? */
439
  if (!cell_is_active_gravity(c, e)) return;
440

Tom Theuns's avatar
Tom Theuns committed
441
442
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
443
    for (int k = 0; k < 8; k++)
444
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
445
  } else {
446

447
448
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
449

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

453
      /* Is this part within the time step? */
454
      if (gpart_is_active(gp, e)) {
455
456
        external_gravity_acceleration(time, potential, constants, gp);
      }
457
    }
458
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
459

460
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
461
462
}

463
464
465
466
467
468
469
470
471
/**
 * @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) {

472
473
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
  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
498
/**
499
500
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
501
502
503
504
505
506
507
 *
 * @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) {

508
  const struct engine *e = r->e;
509
510
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
511
512
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
513
  const struct unit_system *us = e->internal_units;
514
  const struct hydro_props *hydro_props = e->hydro_properties;
515
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
516
  const double time_base = e->time_base;
517
  const integertime_t ti_current = e->ti_current;
518
519
520
  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
521
522
523

  TIMER_TIC;

524
  /* Anything to do here? */
525
  if (!cell_is_active_hydro(c, e)) return;
526

Stefan Arridge's avatar
Stefan Arridge committed
527
528
529
530
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
531
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
532

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

536
537
538
      /* 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
539

540
      if (part_is_active(p, e)) {
541

542
        double dt_cool, dt_therm;
543
544
545
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
546
547
              get_integer_time_begin(ti_current - 1, p->time_bin);

548
549
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
550
551
552
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

553
554
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
555
          dt_therm = get_timestep(p->time_bin, time_base);
556
        }
557

558
        /* Let's cool ! */
559
560
561
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
562
      }
Stefan Arridge's avatar
Stefan Arridge committed
563
564
565
566
567
568
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
569
570
571
572
573
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

574
  struct engine *e = r->e;
575
  const struct cosmology *cosmo = e->cosmology;
576
577
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
578
579
580
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
581
  const int with_cosmology = (e->policy & engine_policy_cosmology);
582
  const int with_feedback = (e->policy & engine_policy_feedback);
583
584
585
  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;
586
  const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
587
588
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
589
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
590
591
592
593

  TIMER_TIC;

  /* Anything to do here? */
594
  if (!cell_is_active_hydro(c, e)) {
595
    star_formation_logger_log_inactive_cell(&c->stars.sfh);
596
597
    return;
  }
Folkert Nobels's avatar
Folkert Nobels committed
598
  star_formation_logger_log_active_cell(&c->stars.sfh);
Matthieu Schaller's avatar
Matthieu Schaller committed
599
600
601
602

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
603
604
605
606
607
608
609
610
      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 */
611
        star_formation_logger_log_progeny_cell(&c->stars.sfh, &cp->stars.sfh);
612
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
613
  } else {
614
615
616
617
618
619
620
621

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

622
      /* Only work on active particles */
623
624
      if (part_is_active(p, e)) {

625
626
        /* 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
627
628
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
629

630
631
632
633
634
635
          /* 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);
636

637
638
639
640
641
642
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

644
645
646
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
647

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

651
652
653
654
655
656
657
658
659
660
          /* 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);

            /* Copy the properties of the gas particle to the star particle */
            star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
                                           with_cosmology);
661
662

            /* Update the Star formation history */
663
            star_formation_logger_log_new_spart(sp, &c->stars.sfh);
664
665
666
667
668
669
670
671
          }

        } 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
672
673
        }      /* Not Star-forming? */
      } else { /* is active? */
674
        /* Check if the particle is not inhibited */
675
        if (!part_is_inhibited(p, e)) {
676
677
          star_formation_logger_log_inactive_part(p, xp, &c->stars.sfh);
        }
678
      }
Folkert Nobels's avatar
Folkert Nobels committed
679
    } /* Loop over particles */
Matthieu Schaller's avatar
Matthieu Schaller committed
680
681
  }

682
683
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
684
685
  if (with_feedback && (c == c->hydro.super) &&
      (current_stars_count != c->stars.count)) {
686
687
    cell_clear_stars_sort_flags(c, /*is_super=*/1);
    runner_do_stars_sort(r, c, 0x1FFF, /*cleanup=*/0, /*timer=*/0);
688
  }
689

Matthieu Schaller's avatar
Matthieu Schaller committed
690
691
692
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
693
694
695
696
697
698
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
699
void runner_do_sort_ascending(struct entry *sort, int N) {
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753

  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
754
        }
755
756
757
758
759
760
761
762
763
764
765
766
      } 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
767
    }
768
769
770
  }
}

771
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
772
773
774
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
775
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
776
 */
777
778
779
780
781
782
783
784
785
#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
786
#else
Loic Hausammann's avatar
Loic Hausammann committed
787
788
789
790
#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
791
#endif
Loic Hausammann's avatar
Loic Hausammann committed
792
793
794

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
795

Pedro Gonnet's avatar
Pedro Gonnet committed
796
797
798
799
800
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
801
 * @param flags Cell flag.
802
803
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
804
805
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
806
 */
Loic Hausammann's avatar
Loic Hausammann committed
807
808
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
809
810

  struct entry *fingers[8];
811
812
813
  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
814
  float buff[8];
815

816
817
  TIMER_TIC;

818
  /* We need to do the local sorts plus whatever was requested further up. */
819
  flags |= c->hydro.do_sort;
820
  if (cleanup) {
821
    c->hydro.sorted = 0;
822
  } else {
823
    flags &= ~c->hydro.sorted;
824
  }
825
  if (flags == 0 && !c->hydro.do_sub_sort) return;
826
827

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

831
832
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
833
  runner_check_sorts_hydro(c, c->hydro.sorted);
834
835

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
836
837
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
838
839
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
840
  }
841

842
843
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
844
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
845
#endif
846

847
848
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
849
850
851
852
853

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

    /* Fill in the gaps within the progeny. */
854
    float dx_max_sort = 0.0f;
855
    float dx_max_sort_old = 0.0f;
856
    for (int k = 0; k < 8; k++) {
857
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
858
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
859
        runner_do_hydro_sort(r, c->progeny[k], flags,
860
861
862
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
863
864
865
        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);
866
      }
867
    }
868
869
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
870
871

    /* Loop over the 13 different sort arrays. */
872
    for (int j = 0; j < 13; j++) {
873
874
875
876
877

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

      /* Init the particle index offsets. */
878
      int off[8];
879
880
      off[0] = 0;
      for (int k = 1; k < 8; k++)
881
        if (c->progeny[k - 1] != NULL)
882
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
883
884
885
886
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
887
      int inds[8];
888
      for (int k = 0; k < 8; k++) {
889
        inds[k] = k;
890
891
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
892
893
894
895
896
897
898
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
899
900
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
901
          if (buff[inds[k]] < buff[inds[i]]) {
902
            int temp_i = inds[i];
903
904
905
906
907
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
908
      struct entry *finger = c->hydro.sort[j];
909
      for (int ind = 0; ind < count; ind++) {
910
911
912
913
914
915
916
917
918
919

        /* 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. */
920
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
921
          int temp_i = inds[k - 1];
922
923
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
924
        }
925

926
927
928
      } /* Merge. */

      /* Add a sentinel. */
929
930
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
931
932

      /* Mark as sorted. */
933
      atomic_or(&c->hydro.sorted, 1 << j);
934
935
936
937
938
939
940
941

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

942
    /* Reset the sort distance */
943
    if (c->hydro.sorted == 0) {
944
945
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
946
        error("Have non-NULL xparts in foreign cell");
947
#endif
948
949
950
951
952
953
954
955

      /* And the individual sort distances if we are a local cell */
      if (xparts != NULL) {
        for (int k = 0; k < count; k++) {
          xparts[k].x_diff_sort[0] = 0.0f;
          xparts[k].x_diff_sort[1] = 0.0f;
          xparts[k].x_diff_sort[2] = 0.0f;
        }
956
      }
957
958
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
959
960
    }

961
    /* Fill the sort array. */
962
    for (int k = 0; k < count; k++) {
963
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
964
      for (int j = 0; j < 13; j++)
965
        if (flags & (1 << j)) {
966
967
968
969
          c->hydro.sort[j][k].i = k;
          c->hydro.sort[j][k].d = px[0] * runner_shift[j][0] +
                                  px[1] * runner_shift[j][1] +
                                  px[2] * runner_shift[j][2];
970