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 "minmax.h"
James Willis's avatar
James Willis committed
59
#include "runner_doiact_vec.h"
60
#include "scheduler.h"
61
#include "sort_part.h"
62
#include "space.h"
63
#include "space_getsid.h"
64
#include "star_formation.h"
65
#include "star_formation_iact.h"
66
#include "stars.h"
67
68
#include "task.h"
#include "timers.h"
69
#include "timestep.h"
70
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
71
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
72

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

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

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

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

102
103
104
105
106
107
108
/* 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

109
/* Import the gravity loop functions. */
110
#include "runner_doiact_grav.h"
111

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

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

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

136
  struct spart *restrict sparts = c->stars.parts;
137
138
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
139
140
141
142
143
144
  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;
145
  int redo = 0, scount = 0;
146
147
148
149

  TIMER_TIC;

  /* Anything to do here? */
Loic Hausammann's avatar
Loic Hausammann committed
150
  if (!cell_is_active_stars(c, e)) return;
151
152
153
154

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

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

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

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

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

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

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

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

205
206
207
        float h_new;
        int has_no_neighbours = 0;

208
        if (sp->density.wcount == 0.f) { /* No neighbours case */
209
210
211
212
213
214

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

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

216
217
218
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
219
          stars_end_density(sp, cosmo);
220
221

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

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

#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

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

246
247
            stars_reset_feedback(sp);

248
249
250
251
252
            /* Ok, we are done with this particle */
            continue;
          }

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

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

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

268
269
270
271
272
273
274
275
276
#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);
277
278
279
280

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

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

          /* Ok, correct then */
287
288
289
290
291
292
293
294
295
296
297
298
299
300

          /* 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;
          }
301
302

          /* If below the absolute maximum, try again */
303
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
304
305
306

            /* Flag for another round of fun */
            sid[redo] = sid[i];
307
308
309
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
310
311
312
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
313
            stars_init_spart(sp);
314
315
316
317

            /* Off we go ! */
            continue;

318
319
320
321
322
323
          } 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) {
324
325

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

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

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
337
338
339
          }
        }

340
        /* We now have a particle whose smoothing length has converged */
341
        stars_reset_feedback(sp);
342

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

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

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

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

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

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

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

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

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

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

403
404
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
405
406
407
    }

    /* Be clean */
408
409
    free(left);
    free(right);
410
    free(sid);
411
    free(h_0);
412
413
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
414
  if (timer) TIMER_TOC(timer_dostars_ghost);
415
416
}

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

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

433
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
434

435
  /* Anything to do here? */
436
  if (!cell_is_active_gravity(c, e)) return;
437

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

444
445
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
446

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

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

457
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
458
459
}

460
461
462
463
464
465
466
467
468
/**
 * @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) {

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

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

  TIMER_TIC;

521
  /* Anything to do here? */
522
  if (!cell_is_active_hydro(c, e)) return;
523

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

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

533
534
535
      /* 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
536

537
      if (part_is_active(p, e)) {
538

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

545
546
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
547
548
549
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

550
551
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
552
          dt_therm = get_timestep(p->time_bin, time_base);
553
        }
554

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
566
567
568
569
570
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

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

  TIMER_TIC;

590
591
592
  /* Initialize the star formation tracers */
  starformation_init_SFH(c->stars.sfh, cosmo, with_cosmology);

Matthieu Schaller's avatar
Matthieu Schaller committed
593
594
595
596
597
598
  /* Anything to do here? */
  if (!cell_is_active_hydro(c, e)) return;

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

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

618
      /* Only work on active particles */
619
620
      if (part_is_active(p, e)) {

621
622
        /* 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
623
624
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
625

626
627
628
629
630
631
          /* 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);
632

633
634
635
636
637
638
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

640
641
642
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
643

644
645
646
647
648
649
650
651
652
653
          /* 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);
654
655
656

            /* Update the Star formation history */
            starformation_update_SFH(sp, c->stars.sfh, cosmo, with_cosmology);
657
658
659
660
661
662
663
664
665
666
667
          }

        } 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
668
669
  }

670
671
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
672
673
  if (with_feedback && (c == c->hydro.super) &&
      (current_stars_count != c->stars.count)) {
674
675
    cell_clear_stars_sort_flags(c, /*is_super=*/1);
    runner_do_stars_sort(r, c, 0x1FFF, /*cleanup=*/0, /*timer=*/0);
676
  }
677

Matthieu Schaller's avatar
Matthieu Schaller committed
678
679
680
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
681
682
683
684
685
686
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
687
void runner_do_sort_ascending(struct entry *sort, int N) {
688
689
690
691
692
693
694
695
696
697
698
699
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

  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
742
        }
743
744
745
746
747
748
749
750
751
752
753
754
      } 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
755
    }
756
757
758
  }
}

759
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
760
761
762
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
763
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
764
 */
765
766
767
768
769
770
771
772
773
#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
774
#else
Loic Hausammann's avatar
Loic Hausammann committed
775
776
777
778
#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
779
#endif
Loic Hausammann's avatar
Loic Hausammann committed
780
781
782

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
783

Pedro Gonnet's avatar
Pedro Gonnet committed
784
785
786
787
788
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
789
 * @param flags Cell flag.
790
791
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
792
793
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
794
 */
Loic Hausammann's avatar
Loic Hausammann committed
795
796
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
797
798

  struct entry *fingers[8];
799
800
801
  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
802
  float buff[8];
803

804
805
  TIMER_TIC;

806
  /* We need to do the local sorts plus whatever was requested further up. */
807
  flags |= c->hydro.do_sort;
808
  if (cleanup) {
809
    c->hydro.sorted = 0;
810
  } else {
811
    flags &= ~c->hydro.sorted;
812
  }
813
  if (flags == 0 && !c->hydro.do_sub_sort) return;
814
815

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

819
820
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
821
  runner_check_sorts_hydro(c, c->hydro.sorted);
822
823

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
824
825
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
826
827
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
828
  }
829

830
831
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
832
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
833
#endif
834

835
836
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
837
838
839
    if ((flags & (1 << j)) && c->hydro.sort[j] == NULL) {
      if ((c->hydro.sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                                     (count + 1))) == NULL)
840
841
        error("Failed to allocate sort memory.");
    }
842
843
844
845
846
847
  }

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

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

    /* Loop over the 13 different sort arrays. */
866
    for (int j = 0; j < 13; j++) {
867
868
869
870
871

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

      /* Init the particle index offsets. */
872
      int off[8];
873
874
      off[0] = 0;
      for (int k = 1; k < 8; k++)
875
        if (c->progeny[k - 1] != NULL)
876
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
877
878
879
880
        else
          off[k] = off[k - 1];

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

      /* Sort the buffer. */
893
894
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
895
          if (buff[inds[k]] < buff[inds[i]]) {
896
            int temp_i = inds[i];
897
898
899
900
901
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
902
      struct entry *finger = c->hydro.sort[j];
903
      for (int ind = 0; ind < count; ind++) {
904
905
906
907
908
909
910
911
912
913

        /* 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. */
914
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
915
          int temp_i = inds[k - 1];
916
917
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
918
        }
919

920
921
922
      } /* Merge. */

      /* Add a sentinel. */
923
924
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
925
926

      /* Mark as sorted. */
927
      atomic_or(&c->hydro.sorted, 1 << j);
928
929
930
931
932
933
934
935

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

936
    /* Reset the sort distance */
937
    if (c->hydro.sorted == 0) {
938
939
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
940
        error("Have non-NULL xparts in foreign cell");
941
#endif
942
943
944
945
946
947
948
949

      /* 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;
        }
950
      }
951
952
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
953
954
    }

955
    /* Fill the sort array. */
956
    for (int k = 0; k < count; k++) {
957
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
958
      for (int j = 0; j < 13; j++)
959
        if (flags & (1 << j)) {
960
961
962
963
          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];
964
        }
965
    }
966
967

    /* Add the sentinel and sort. */
968
    for (int j = 0; j < 13; j++)
969
      if (flags & (1 << j)) {
970
971
972
973
        c->hydro.sort[j][count].d = FLT_MAX;
        c->hydro.sort[j][count].i = 0;
        runner_do_sort_ascending(c->hydro.sort[j], count);
        atomic_or(&c->hydro.sorted, 1 << j);
974
975
976
      }
  }

977
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
978
  /* Verify the sorting. */
979
  for (int j = 0; j < 13; j++) {
980
    if (!(flags & (1 << j))) continue;
981
    struct entry *finger = c->hydro.sort[j];
982
    for (int k = 1; k < count; k++)