runner.c 111 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 "star_formation_logger.h"
67
#include "stars.h"
68
69
#include "task.h"
#include "timers.h"
70
#include "timestep.h"
71
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
72
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
73

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

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

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

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

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

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

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

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

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

206
207
208
        float h_new;
        int has_no_neighbours = 0;

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

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

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

217
218
219
        } else {

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

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

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

#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

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

247
248
            stars_reset_feedback(sp);

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

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

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

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

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

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

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

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

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

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

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

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

            /* Off we go ! */
            continue;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

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

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

  TIMER_TIC;

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

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

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

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

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

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

635
636
637
638
639
640
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

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

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

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

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

        } 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
670
671
        }      /* Not Star-forming? */
      } else { /* is active? */
672
        star_formation_logger_log_inactive_part(p, xp, &c->stars.sfh);
673
      }
Folkert Nobels's avatar
Folkert Nobels committed
674
    } /* Loop over particles */
Matthieu Schaller's avatar
Matthieu Schaller committed
675
676
  }

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

Matthieu Schaller's avatar
Matthieu Schaller committed
685
686
687
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
688
689
690
691
692
693
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
694
void runner_do_sort_ascending(struct entry *sort, int N) {
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
742
743
744
745
746
747
748

  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
749
        }
750
751
752
753
754
755
756
757
758
759
760
761
      } 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
762
    }
763
764
765
  }
}

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

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
790

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

  struct entry *fingers[8];
806
807
808
  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
809
  float buff[8];
810

811
812
  TIMER_TIC;

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

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

826
827
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
828
  runner_check_sorts_hydro(c, c->hydro.sorted);
829
830

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

837
838
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
839
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
840
#endif
841

842
843
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
844
845
846
    if ((flags & (1 << j)) && c->hydro.sort[j] == NULL) {
      if ((c->hydro.sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                                     (count + 1))) == NULL)
847
848
        error("Failed to allocate sort memory.");
    }
849
850
851
852
853
854
  }

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

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

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

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

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

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

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

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

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

927
928
929
      } /* Merge. */

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

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

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

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

      /* 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;
        }
957
      }
958
959
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
960
961
    }

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