runner.c 121 KB
Newer Older
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
5
6
7
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
 *               2016 John A. Regan (john.a.regan@durham.ac.uk)
 *                    Tom Theuns (tom.theuns@durham.ac.uk)
8
 *
9
10
11
12
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
13
 *
14
15
16
17
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
18
 *
19
20
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
21
 *
22
 ******************************************************************************/
Pedro Gonnet's avatar
Pedro Gonnet committed
23

Pedro Gonnet's avatar
Pedro Gonnet committed
24
25
/* Config parameters. */
#include "../config.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
26
27
28
29

/* Some standard headers. */
#include <float.h>
#include <limits.h>
30
#include <stdlib.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
31

32
33
/* MPI headers. */
#ifdef WITH_MPI
34
#include <mpi.h>
35
36
#endif

37
38
39
/* This object's header. */
#include "runner.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
40
/* Local headers. */
41
#include "active.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
42
#include "approx_math.h"
43
#include "atomic.h"
44
#include "cell.h"
45
#include "chemistry.h"
46
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
47
#include "cooling.h"
48
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
49
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
50
#include "engine.h"
51
#include "entropy_floor.h"
52
#include "error.h"
53
#include "feedback.h"
54
55
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
56
#include "hydro_properties.h"
57
#include "kick.h"
Loikki's avatar
Loikki committed
58
#include "logger.h"
59
#include "memuse.h"
60
#include "minmax.h"
James Willis's avatar
James Willis committed
61
#include "runner_doiact_vec.h"
62
#include "scheduler.h"
63
#include "sort_part.h"
64
#include "space.h"
65
#include "space_getsid.h"
66
#include "star_formation.h"
67
#include "star_formation_iact.h"
68
#include "star_formation_logger.h"
69
#include "stars.h"
70
71
#include "task.h"
#include "timers.h"
72
#include "timestep.h"
73
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
74
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
75

76
77
78
79
#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
80
#define TASK_LOOP_FEEDBACK 4
81

82
/* Import the density loop functions. */
83
#define FUNCTION density
84
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
85
#include "runner_doiact.h"
86
87
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
88

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

98
/* Import the force loop functions. */
99
#define FUNCTION force
100
#define FUNCTION_TASK_LOOP TASK_LOOP_FORCE
101
#include "runner_doiact.h"
102
103
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
104

105
106
107
108
109
110
111
/* 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

112
/* Import the gravity loop functions. */
113
#include "runner_doiact_grav.h"
114

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

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

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

139
  struct spart *restrict sparts = c->stars.parts;
140
  const struct engine *e = r->e;
141
142
  const struct unit_system *us = e->internal_units;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
143
  const struct cosmology *cosmo = e->cosmology;
144
  const struct feedback_props *feedback_props = e->feedback_props;
145
146
147
148
149
150
  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;
151
  int redo = 0, scount = 0;
152

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

156
157
158
  TIMER_TIC;

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

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

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

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

    /* While there are particles that need to be updated... */
198
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
199
200
201
202
203
204
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
205
      for (int i = 0; i < scount; i++) {
206
207
208
209
210
211

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
212
213
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
214
215
216
#endif

        /* Get some useful values */
217
        const float h_init = h_0[i];
218
219
220
        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);
221

222
223
224
        float h_new;
        int has_no_neighbours = 0;

225
        if (sp->density.wcount == 0.f) { /* No neighbours case */
226
227
228
229
230
231

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

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

233
234
235
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
236
          stars_end_density(sp, cosmo);
237
238

          /* Compute one step of the Newton-Raphson scheme */
239
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
240
          const float n_target = stars_eta_dim;
241
242
          const float f = n_sum - n_target;
          const float f_prime =
243
244
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
245

246
          /* Improve the bisection bounds */
247
248
249
250
          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);
251
252
253
254
255
256
257

#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

258
          /* Skip if h is already h_max and we don't have enough neighbours */
259
260
261
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
262

263
            stars_reset_feedback(sp);
264
            feedback_reset_feedback(sp, feedback_props);
265

266
267
268
269
270
            /* Ok, we are done with this particle */
            continue;
          }

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

272
273
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
274
275
276
277
278
279
280
281
282
283
284
285

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

286
287
288
          /* 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
            feedback_init_spart(sp);
327
328
329
330

            /* Off we go ! */
            continue;

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

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

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

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

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

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

358
        stars_reset_feedback(sp);
359

360
        /* Only do feedback if stars have a reasonable birth time */
361
        if (feedback_do_feedback(sp)) {
362

363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
          const integertime_t ti_step = get_integer_timestep(sp->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(e->ti_current - 1, sp->time_bin);

          /* Get particle time-step */
          double dt;
          if (with_cosmology) {
            dt = cosmology_get_delta_time(e->cosmology, ti_begin,
                                          ti_begin + ti_step);
          } else {
            dt = get_timestep(sp->time_bin, e->time_base);
          }

          /* Calculate age of the star at current time */
          double star_age_end_of_step;
378
          if (with_cosmology) {
379
            star_age_end_of_step = cosmology_get_delta_time_from_scale_factors(
380
381
                cosmo, sp->birth_scale_factor, (float)cosmo->a);
          } else {
382
            star_age_end_of_step = (float)e->time - sp->birth_time;
383
384
          }

385
386
387
388
389
390
391
392
393
394
          /* Has this star been around for a while ? */
          if (star_age_end_of_step > 0.) {

            /* Age of the star at the start of the step */
            const double star_age_beg_of_step =
                max(star_age_end_of_step - dt, 0.);

            /* Compute the stellar evolution  */
            feedback_evolve_spart(sp, feedback_props, cosmo, us,
                                  star_age_beg_of_step, dt);
395
396
397
398
          } else {

            /* Reset the feedback fields of the star particle */
            feedback_reset_feedback(sp, feedback_props);
399
          }
400
401
402
403
        } else {

          /* Reset the feedback fields of the star particle */
          feedback_reset_feedback(sp, feedback_props);
404
        }
405
406
407
408
409
410
      }

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

      /* Re-set the counter for the next loop (potentially). */
411
412
      scount = redo;
      if (scount > 0) {
413
414
415
416
417

        /* 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
418
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
419
420
421
422
423
424
425
426

#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)
427
428
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
429
430
431
432
433
434

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

              /* Left or right? */
              if (l->t->ci == finger)
435
436
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
437
              else
438
439
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
440
441
442
443
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
444
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
445
                                                NULL, 1);
446
447
448
449
450
451

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

              /* Left or right? */
              if (l->t->ci == finger)
452
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
453
                                                  scount, l->t->cj, 1);
454
              else
455
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
456
                                                  scount, l->t->ci, 1);
457
458
459
460
461
462
            }
          }
        }
      }
    }

463
464
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
465
466
467
    }

    /* Be clean */
468
469
    free(left);
    free(right);
470
    free(sid);
471
    free(h_0);
472
473
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
474
475
  /* Update h_max */
  c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
476

477
  /* The ghost may not always be at the top level.
478
   * Therefore we need to update h_max between the super- and top-levels */
479
  if (c->stars.ghost) {
480
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
481
      atomic_max_d(&tmp->stars.h_max, h_max);
482
483
484
    }
  }

485
  if (timer) TIMER_TOC(timer_do_stars_ghost);
486
487
}

Tom Theuns's avatar
Tom Theuns committed
488
489
490
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
491
492
493
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
494
 */
495
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
496

497
498
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
499
500
501
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
502
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
503

504
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
505

506
  /* Anything to do here? */
507
  if (!cell_is_active_gravity(c, e)) return;
508

Tom Theuns's avatar
Tom Theuns committed
509
510
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
511
    for (int k = 0; k < 8; k++)
512
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
513
  } else {
514

515
516
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
517

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

521
      /* Is this part within the time step? */
522
      if (gpart_is_active(gp, e)) {
523
524
        external_gravity_acceleration(time, potential, constants, gp);
      }
525
    }
526
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
527

528
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
529
530
}

531
532
533
534
535
536
537
538
539
/**
 * @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) {

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

576
  const struct engine *e = r->e;
577
578
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
579
580
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
581
  const struct unit_system *us = e->internal_units;
582
  const struct hydro_props *hydro_props = e->hydro_properties;
583
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
584
  const double time_base = e->time_base;
585
  const integertime_t ti_current = e->ti_current;
586
587
588
  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
589
590
591

  TIMER_TIC;

592
  /* Anything to do here? */
593
  if (!cell_is_active_hydro(c, e)) return;
594

Stefan Arridge's avatar
Stefan Arridge committed
595
596
597
598
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
599
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
600

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

604
605
606
      /* 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
607

608
      if (part_is_active(p, e)) {
609

610
        double dt_cool, dt_therm;
611
612
613
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
614
615
              get_integer_time_begin(ti_current - 1, p->time_bin);

616
617
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
618
619
620
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

621
622
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
623
          dt_therm = get_timestep(p->time_bin, time_base);
624
        }
625

626
        /* Let's cool ! */
627
628
629
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
630
      }
Stefan Arridge's avatar
Stefan Arridge committed
631
632
633
634
635
636
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
637
638
639
640
641
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

642
  struct engine *e = r->e;
643
  const struct cosmology *cosmo = e->cosmology;
644
645
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
646
647
648
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
649
  const int with_cosmology = (e->policy & engine_policy_cosmology);
650
  const int with_feedback = (e->policy & engine_policy_feedback);
651
652
653
  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;
654
  const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
655
656
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
657
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
658
659
660

  TIMER_TIC;

661
662
663
664
665
#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
666
  /* Anything to do here? */
667
  if (c->hydro.count == 0 || !cell_is_active_hydro(c, e)) {
668
    star_formation_logger_log_inactive_cell(&c->stars.sfh);
669
670
    return;
  }
671
672

  /* Reset the SFR */
673
  star_formation_logger_init(&c->stars.sfh);
Matthieu Schaller's avatar
Matthieu Schaller committed
674
675
676
677

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
678
679
680
681
682
683
684
685
      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 */
686
        star_formation_logger_add(&c->stars.sfh, &cp->stars.sfh);
687
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
688
  } else {
689
690
691
692
693
694
695
696

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

697
      /* Only work on active particles */
698
699
      if (part_is_active(p, e)) {

700
701
        /* 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
702
703
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
704

705
706
707
708
709
710
          /* 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);
711

712
713
714
715
716
717
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

719
720
721
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
722

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

726
727
728
729
730
731
732
          /* 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);

733
734
735
736
737
            /* 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,
738
739
                                             with_cosmology, phys_const, hydro_props,
                                             us, cooling);
740
741
742

              /* Update the Star formation history */
              star_formation_logger_log_new_spart(sp, &c->stars.sfh);
743
            }
744
745
746
747
748
749
750
751
752
          }

        } 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? */
753

Folkert Nobels's avatar
Folkert Nobels committed
754
      } else { /* is active? */
755

756
        /* Check if the particle is not inhibited */
757
        if (!part_is_inhibited(p, e)) {
758
759
          star_formation_logger_log_inactive_part(p, xp, &c->stars.sfh);
        }
760
      }
Folkert Nobels's avatar
Folkert Nobels committed
761
    } /* Loop over particles */
Matthieu Schaller's avatar
Matthieu Schaller committed
762
763
  }

764
765
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
766
  if (with_feedback && (c == c->top) &&
767
      (current_stars_count != c->stars.count)) {
768
769

    cell_clear_stars_sort_flags(c);
770
    runner_do_all_stars_sort(r, c);
771
  }
772

Matthieu Schaller's avatar
Matthieu Schaller committed
773
774
775
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
776
777
778
779
780
781
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
782
void runner_do_sort_ascending(struct entry *sort, int N) {
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836

  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
837
        }
838
839
840
841
842
843
844
845
846
847
848
849
      } 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
850
    }
851
852
853
  }
}

854
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
855
856
857
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
858
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
859
 */
860
861
862
863
864
865
866
867
868
#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
869
#else
Loic Hausammann's avatar
Loic Hausammann committed
870
871
872
873
#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
874
#endif
Loic Hausammann's avatar
Loic Hausammann committed
875
876
877

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
878

Pedro Gonnet's avatar
Pedro Gonnet committed
879
880
881
882
883
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
884
 * @param flags Cell flag.
885
886
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
887
888
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
889
 */
Loic Hausammann's avatar
Loic Hausammann committed
890
891
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
892
893

  struct entry *fingers[8];
894
895
896
  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
897
  float buff[8];
898

899
900
  TIMER_TIC;

901
902
903
904
#ifdef SWIFT_DEBUG_CHECKS
  if (c->hydro.super == NULL) error("Task called above the super level!!!");
#endif

905
  /* We need to do the local sorts plus whatever was requested further up. */
906
  flags |= c->hydro.do_sort;
907
  if (cleanup) {
908
    c->hydro.sorted = 0;
909
  } else {
910
    flags &= ~c->hydro.sorted;
911
  }
912
  if (flags == 0 && !c->hydro.do_sub_sort) return;
913
914

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

918
919
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
920
  runner_check_sorts_hydro(c, c->hydro.sorted);
921
922

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
923
924
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
925
926
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
927
  }
928

929
930
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
931
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
932
#endif
933

934
935
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
936
937
938
939
940

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

    /* Fill in the gaps within the progeny. */
941
    float dx_max_sort = 0.0f;
942
    float dx_max_sort_old = 0.0f;
943
    for (int k = 0; k < 8; k++) {
944
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
945
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
946
        runner_do_hydro_sort(r, c->progeny[k], flags,
947
948
949
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
950
951
952
        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);
953
      }
954
    }
955
956
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
957
958

    /* Loop over the 13 different sort arrays. */
959
    for (int j = 0; j < 13; j++) {