runner.c 120 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"
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
      if (spart_is_active(&sparts[k], e)) {
189
        sid[scount] = k;
190
191
192
        h_0[scount] = sparts[k].h;
        left[scount] = 0.f;
        right[scount] = stars_h_max;
193
        ++scount;
194
195
196
      }

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

200
201
202
203
      /* Reset the redo-count. */
      redo = 0;

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

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

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

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

221
222
223
        float h_new;
        int has_no_neighbours = 0;

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

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

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

232
233
234
        } else {

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

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

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

#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

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

262
            stars_reset_feedback(sp);
263
264
265
266
267
268

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

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

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

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

284
285
286
          /* 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);
287
288
289
290

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

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

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

          /* 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;
          }
311
312

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

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

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

            /* Off we go ! */
            continue;

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

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

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

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

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

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

356
        stars_reset_feedback(sp);
357

358
        /* Only do feedback if stars have a reasonable birth time */
359
360
        if (sp->birth_time != -1.) {

361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
          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;
376
          if (with_cosmology) {
377
            star_age_end_of_step = cosmology_get_delta_time_from_scale_factors(
378
379
                cosmo, sp->birth_scale_factor, (float)cosmo->a);
          } else {
380
            star_age_end_of_step = e->time - sp->birth_time;
381
382
          }

383
384
385
386
387
388
389
390
391
392
393
394
395
396
          /* Reset the feedback fields of the star particle */
          feedback_prepare_spart(sp, feedback_props);

          /* 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);
          }
397
        }
398
399
400
401
402
403
      }

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

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

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

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

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

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

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

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

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

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

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

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
478
  if (timer) TIMER_TOC(timer_dostars_ghost);
479
480
}

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

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

497
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
498

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

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

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

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

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

521
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
522
523
}

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

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

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

  TIMER_TIC;

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

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

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

597
598
599
      /* 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
600

601
      if (part_is_active(p, e)) {
602

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

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

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

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

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

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

  TIMER_TIC;

654
655
656
657
658
#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
659
  /* Anything to do here? */
660
  if (c->hydro.count == 0) return;
661
  if (!cell_is_active_hydro(c, e)) {
662
    star_formation_logger_log_inactive_cell(&c->stars.sfh);
663
664
    return;
  }
665
  star_formation_logger_init(&c->stars.sfh);
Matthieu Schaller's avatar
Matthieu Schaller committed
666
667
668
669

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
670
671
672
673
674
675
676
677
      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 */
678
        star_formation_logger_add(&c->stars.sfh, &cp->stars.sfh);
679
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
680
  } else {
681
682
683
684
685
686
687
688

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

689
      /* Only work on active particles */
690
691
      if (part_is_active(p, e)) {

692
693
        /* 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
694
695
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
696

697
698
699
700
701
702
          /* 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);
703

704
705
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
706

707
708
709
          } else {
            dt_star = get_timestep(p->time_bin, time_base);
          }
710

711
712
713
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
714

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

718
719
720
721
722
723
724
          /* 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);

725
726
727
            /* Did we get a star? (Or did we run out of spare ones?) */
            if (sp != NULL) {

728
729
              message("Formed a star ID=%lld", sp->id);

730
731
732
              /* Copy the properties of the gas particle to the star particle */
              star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
                                             with_cosmology);
733
734
735

              /* Update the Star formation history */
              star_formation_logger_log_new_spart(sp, &c->stars.sfh);
736
            }
737
738
739
740
741
742
743
744
          }

        } 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
745
746
        }      /* Not Star-forming? */
      } else { /* is active? */
747
        /* Check if the particle is not inhibited */
748
        if (!part_is_inhibited(p, e)) {
749
750
          star_formation_logger_log_inactive_part(p, xp, &c->stars.sfh);
        }
751
      }
Folkert Nobels's avatar
Folkert Nobels committed
752
    } /* Loop over particles */
Matthieu Schaller's avatar
Matthieu Schaller committed
753
754
  }

755
756
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
757
  if (with_feedback && (c == c->top) &&
758
      (current_stars_count != c->stars.count)) {
759
760

    cell_clear_stars_sort_flags(c);
761
    runner_do_all_stars_sort(r, c);
762
  }
763

Matthieu Schaller's avatar
Matthieu Schaller committed
764
765
766
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
767
768
769
770
771
772
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
773
void runner_do_sort_ascending(struct entry *sort, int N) {
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
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

  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
828
        }
829
830
831
832
833
834
835
836
837
838
839
840
      } 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
841
    }
842
843
844
  }
}

845
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
846
847
848
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
849
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
850
 */
851
852
853
854
855
856
857
858
859
#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
860
#else
Loic Hausammann's avatar
Loic Hausammann committed
861
862
863
864
#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
865
#endif
Loic Hausammann's avatar
Loic Hausammann committed
866
867
868

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
869

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

  struct entry *fingers[8];
885
886
887
  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
888
  float buff[8];
889

890
891
  TIMER_TIC;

892
893
894
895
#ifdef SWIFT_DEBUG_CHECKS
  if (c->hydro.super == NULL) error("Task called above the super level!!!");
#endif

896
  /* We need to do the local sorts plus whatever was requested further up. */
897
  flags |= c->hydro.do_sort;
898
  if (cleanup) {
899
    c->hydro.sorted = 0;
900
  } else {
901
    flags &= ~c->hydro.sorted;
902
  }
903
  if (flags == 0 && !c->hydro.do_sub_sort) return;
904
905

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

909
910
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
911
  runner_check_sorts_hydro(c, c->hydro.sorted);
912
913

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
914
915
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
916
917
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
918
  }
919

920
921
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
922
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
923
#endif
924

925
926
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
927
928
929
930
931

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

    /* Fill in the gaps within the progeny. */
932
    float dx_max_sort = 0.0f;
933
    float dx_max_sort_old = 0.0f;
934
    for (int k = 0; k < 8; k++) {
935
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
936
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
937
        runner_do_hydro_sort(r, c->progeny[k], flags,
938
939
940
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
941
942
943
        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);
944
      }
945
    }
946
947
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
948
949

    /* Loop over the 13 different sort arrays. */
950
    for (int j = 0; j < 13; j++) {
951
952
953
954
955

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

      /* Init the particle index offsets. */
956
      int off[8];
957
958
      off[0] = 0;
      for (int k = 1; k < 8; k++)
959
        if (c->progeny[k - 1] != NULL)
960
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;