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

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

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

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

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

Pedro Gonnet's avatar
Pedro Gonnet committed
40
/* Local headers. */
41
#include "active.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
42
#include "approx_math.h"
43
#include "atomic.h"
44
#include "cell.h"
45
#include "chemistry.h"
46
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
47
#include "cooling.h"
48
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
49
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
50
#include "engine.h"
51
#include "entropy_floor.h"
52
#include "error.h"
53
54
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
55
#include "hydro_properties.h"
56
#include "kick.h"
Loikki's avatar
Loikki committed
57
#include "logger.h"
58
#include "memuse.h"
59
#include "minmax.h"
James Willis's avatar
James Willis committed
60
#include "runner_doiact_vec.h"
61
#include "scheduler.h"
62
#include "sort_part.h"
63
#include "space.h"
64
#include "space_getsid.h"
65
#include "star_formation.h"
66
#include "star_formation_iact.h"
67
#include "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

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

151
152
153
  TIMER_TIC;

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
163
164
        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->stars.h_max);
165
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
166
    }
167
168
169
170
  } else {

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

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

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

      /* Loop over the remaining active parts in this cell. */
199
      for (int i = 0; i < scount; i++) {
200
201
202
203
204
205

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
206
207
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
208
209
210
#endif

        /* Get some useful values */
211
        const float h_init = h_0[i];
212
213
214
        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);
215

216
217
218
        float h_new;
        int has_no_neighbours = 0;

219
        if (sp->density.wcount == 0.f) { /* No neighbours case */
220
221
222
223
224
225

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

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

227
228
229
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
230
          stars_end_density(sp, cosmo);
231
232

          /* Compute one step of the Newton-Raphson scheme */
233
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
234
          const float n_target = stars_eta_dim;
235
236
          const float f = n_sum - n_target;
          const float f_prime =
237
238
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
239

240
          /* Improve the bisection bounds */
241
242
243
244
          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);
245
246
247
248
249
250
251

#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

252
          /* Skip if h is already h_max and we don't have enough neighbours */
253
254
255
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
256

257
258
            stars_reset_feedback(sp);

259
260
261
262
263
            /* Ok, we are done with this particle */
            continue;
          }

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

265
266
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
267
268
269
270
271
272
273
274
275
276
277
278

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

279
280
281
282
283
284
285
286
287
#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);
288
289
290
291

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

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

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

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

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

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

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
324
            stars_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

Loic Hausammann's avatar
Loic Hausammann committed
358
        /* Compute the stellar evolution  */
359
        stars_evolve_spart(sp, e->stars_properties, cosmo);
360
361
362
363
364
365
      }

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

      /* Re-set the counter for the next loop (potentially). */
366
367
      scount = redo;
      if (scount > 0) {
368
369
370
371
372

        /* 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
373
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
374
375
376
377
378
379
380
381

#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)
382
383
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
384
385
386
387
388
389

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

              /* Left or right? */
              if (l->t->ci == finger)
390
391
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
392
              else
393
394
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
395
396
397
398
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
399
400
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
401
402
403
404
405
406

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

              /* Left or right? */
              if (l->t->ci == finger)
407
408
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
409
              else
410
411
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
412
413
414
415
416
417
            }
          }
        }
      }
    }

418
419
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
420
421
422
    }

    /* Be clean */
423
424
    free(left);
    free(right);
425
    free(sid);
426
    free(h_0);
427
428
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
429
430
  /* Update h_max */
  c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
431

432
  /* The ghost may not always be at the top level.
433
   * Therefore we need to update h_max between the super- and top-levels */
434
  if (c->stars.ghost) {
435
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
436
      atomic_max_f(&tmp->stars.h_max, h_max);
437
438
439
    }
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
440
  if (timer) TIMER_TOC(timer_dostars_ghost);
441
442
}

Tom Theuns's avatar
Tom Theuns committed
443
444
445
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
446
447
448
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
449
 */
450
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
451

452
453
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
454
455
456
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
457
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
458

459
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
460

461
  /* Anything to do here? */
462
  if (!cell_is_active_gravity(c, e)) return;
463

Tom Theuns's avatar
Tom Theuns committed
464
465
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
466
    for (int k = 0; k < 8; k++)
467
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
468
  } else {
469

470
471
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
472

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

476
      /* Is this part within the time step? */
477
      if (gpart_is_active(gp, e)) {
478
479
        external_gravity_acceleration(time, potential, constants, gp);
      }
480
    }
481
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
482

483
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
484
485
}

486
487
488
489
490
491
492
493
494
/**
 * @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) {

495
496
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
  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
521
/**
522
523
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
524
525
526
527
528
529
530
 *
 * @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) {

531
  const struct engine *e = r->e;
532
533
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
534
535
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
536
  const struct unit_system *us = e->internal_units;
537
  const struct hydro_props *hydro_props = e->hydro_properties;
538
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
539
  const double time_base = e->time_base;
540
  const integertime_t ti_current = e->ti_current;
541
542
543
  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
544
545
546

  TIMER_TIC;

547
  /* Anything to do here? */
548
  if (!cell_is_active_hydro(c, e)) return;
549

Stefan Arridge's avatar
Stefan Arridge committed
550
551
552
553
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
554
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
555

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

559
560
561
      /* 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
562

563
      if (part_is_active(p, e)) {
564

565
        double dt_cool, dt_therm;
566
567
568
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
569
570
              get_integer_time_begin(ti_current - 1, p->time_bin);

571
572
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
573
574
575
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

576
577
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
578
          dt_therm = get_timestep(p->time_bin, time_base);
579
        }
580

581
        /* Let's cool ! */
582
583
584
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
585
      }
Stefan Arridge's avatar
Stefan Arridge committed
586
587
588
589
590
591
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
592
593
594
595
596
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

597
  struct engine *e = r->e;
598
  const struct cosmology *cosmo = e->cosmology;
599
600
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
601
602
603
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
604
  const int with_cosmology = (e->policy & engine_policy_cosmology);
605
  const int with_feedback = (e->policy & engine_policy_feedback);
606
607
608
  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;
609
  const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
610
611
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
612
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
613
614
615

  TIMER_TIC;

616
617
618
619
620
#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
621
  /* Anything to do here? */
622
  if (c->hydro.count == 0) return;
Matthieu Schaller's avatar
Matthieu Schaller committed
623
624
625
626
627
628
629
  if (!cell_is_active_hydro(c, e)) return;

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_star_formation(r, c->progeny[k], 0);
  } else {
630
631
632
633
634
635
636
637

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

638
      /* Only work on active particles */
639
640
      if (part_is_active(p, e)) {

641
642
        /* 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
643
644
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
645

646
647
648
649
650
651
          /* 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);
652

653
654
655
656
657
658
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

660
661
662
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
663

664
665
666
667
668
669
670
          /* 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);

671
672
673
674
675
676
677
            /* 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,
                                             with_cosmology);
            }
678
679
680
681
682
683
684
685
686
687
688
          }

        } else { /* Are we not star-forming? */

          /* Update the particle to flag it as not star-forming */
          star_formation_update_part_not_SFR(p, xp, e, sf_props,
                                             with_cosmology);

        } /* Not Star-forming? */
      }   /* is active? */
    }     /* Loop over particles */
Matthieu Schaller's avatar
Matthieu Schaller committed
689
690
  }

691
692
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
693
  if (with_feedback && (c == c->top) &&
694
      (current_stars_count != c->stars.count)) {
695
696

    cell_clear_stars_sort_flags(c);
697
    runner_do_all_stars_sort(r, c);
698
  }
699

Matthieu Schaller's avatar
Matthieu Schaller committed
700
701
702
  if (timer) TIMER_TOC(timer_do_star_formation);
}

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

  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
764
        }
765
766
767
768
769
770
771
772
773
774
775
776
      } 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
777
    }
778
779
780
  }
}

781
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
782
783
784
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
785
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
786
 */
787
788
789
790
791
792
793
794
795
#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
796
#else
Loic Hausammann's avatar
Loic Hausammann committed
797
798
799
800
#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
801
#endif
Loic Hausammann's avatar
Loic Hausammann committed
802
803
804

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
805

Pedro Gonnet's avatar
Pedro Gonnet committed
806
807
808
809
810
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
811
 * @param flags Cell flag.
812
813
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
814
815
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
816
 */
Loic Hausammann's avatar
Loic Hausammann committed
817
818
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
819
820

  struct entry *fingers[8];
821
822
823
  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
824
  float buff[8];
825

826
827
  TIMER_TIC;

828
829
830
831
#ifdef SWIFT_DEBUG_CHECKS
  if (c->hydro.super == NULL) error("Task called above the super level!!!");
#endif

832
  /* We need to do the local sorts plus whatever was requested further up. */
833
  flags |= c->hydro.do_sort;
834
  if (cleanup) {
835
    c->hydro.sorted = 0;
836
  } else {
837
    flags &= ~c->hydro.sorted;
838
  }
839
  if (flags == 0 && !c->hydro.do_sub_sort) return;
840
841

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

845
846
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
847
  runner_check_sorts_hydro(c, c->hydro.sorted);
848
849

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
850
851
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
852
853
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
854
  }
855

856
857
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
858
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
859
#endif
860

861
862
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
863
864
865
866
867

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

    /* Fill in the gaps within the progeny. */
868
    float dx_max_sort = 0.0f;
869
    float dx_max_sort_old = 0.0f;
870
    for (int k = 0; k < 8; k++) {
871
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
872
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
873
        runner_do_hydro_sort(r, c->progeny[k], flags,
874
875
876
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
877
878
879
        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);
880
      }
881
    }
882
883
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
884
885

    /* Loop over the 13 different sort arrays. */
886
    for (int j = 0; j < 13; j++) {
887
888
889
890
891

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

      /* Init the particle index offsets. */
892
      int off[8];
893
894
      off[0] = 0;
      for (int k = 1; k < 8; k++)
895
        if (c->progeny[k - 1] != NULL)
896
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
897
898
899
900
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
901
      int inds[8];
902
      for (int k = 0; k < 8; k++) {
903
        inds[k] = k;
904
905
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
906
907
908
909
910
911
912
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
913
914
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
915
          if (buff[inds[k]] < buff[inds[i]]) {
916
            int temp_i = inds[i];
917
918
919
920
921
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
922
      struct entry *finger = c->hydro.sort[j];
923
      for (int ind = 0; ind < count; ind++) {
924
925
926
927
928
929
930
931
932
933

        /* 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. */
934
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
935
          int temp_i = inds[k - 1];
936
937
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
938
        }
939

940
941
942
      } /* Merge. */

      /* Add a sentinel. */
943
944
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
945
946

      /* Mark as sorted. */
947
      atomic_or(&c->hydro.sorted, 1 << j);
948
949
950
951
952
953
954
955

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

956
    /* Reset the sort distance */
957
    if (c->hydro.sorted == 0) {
958
959
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
960
        error("Have non-NULL xparts in foreign cell");
961
#endif
962
963
964
965
966
967
968
969

      /* 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;
        }
970
      }
971
972
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
973
974
    }

975
    /* Fill the sort array. */