runner.c 161 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 "black_holes.h"
45
#include "black_holes_properties.h"
46
#include "cell.h"
47
#include "chemistry.h"
48
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
49
#include "cooling.h"
50
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
51
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
52
#include "engine.h"
53
#include "entropy_floor.h"
54
#include "error.h"
55
#include "feedback.h"
56
57
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
58
#include "hydro_properties.h"
59
#include "kick.h"
Loikki's avatar
Loikki committed
60
#include "logger.h"
61
#include "memuse.h"
62
#include "minmax.h"
James Willis's avatar
James Willis committed
63
#include "runner_doiact_vec.h"
64
#include "scheduler.h"
65
#include "sort_part.h"
66
#include "space.h"
67
#include "space_getsid.h"
68
#include "star_formation.h"
69
#include "star_formation_iact.h"
70
#include "star_formation_logger.h"
71
#include "stars.h"
72
73
#include "task.h"
#include "timers.h"
74
#include "timestep.h"
75
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
76
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
77

78
/* Unique identifier of loop types */
79
80
81
82
#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
83
#define TASK_LOOP_FEEDBACK 4
84
#define TASK_LOOP_SWALLOW 5
85

86
/* Import the density loop functions. */
87
#define FUNCTION density
88
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
89
#include "runner_doiact.h"
90
91
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
92

93
/* Import the gradient loop functions (if required). */
94
95
#ifdef EXTRA_HYDRO_LOOP
#define FUNCTION gradient
96
#define FUNCTION_TASK_LOOP TASK_LOOP_GRADIENT
97
#include "runner_doiact.h"
98
99
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
100
101
#endif

102
/* Import the force loop functions. */
103
#define FUNCTION force
104
#define FUNCTION_TASK_LOOP TASK_LOOP_FORCE
105
#include "runner_doiact.h"
106
107
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
108

109
110
111
112
113
114
115
/* 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

116
/* Import the gravity loop functions. */
117
#include "runner_doiact_grav.h"
118

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

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

133
134
135
136
137
138
139
/* Import the black hole density loop functions. */
#define FUNCTION density
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
#include "runner_doiact_black_holes.h"
#undef FUNCTION_TASK_LOOP
#undef FUNCTION

140
141
142
143
144
145
146
/* Import the black hole feedback loop functions. */
#define FUNCTION swallow
#define FUNCTION_TASK_LOOP TASK_LOOP_SWALLOW
#include "runner_doiact_black_holes.h"
#undef FUNCTION_TASK_LOOP
#undef FUNCTION

147
148
149
150
151
152
153
/* Import the black hole feedback loop functions. */
#define FUNCTION feedback
#define FUNCTION_TASK_LOOP TASK_LOOP_FEEDBACK
#include "runner_doiact_black_holes.h"
#undef FUNCTION_TASK_LOOP
#undef FUNCTION

154
155
156
157
158
159
160
161
/**
 * @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
162
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
163

164
  struct spart *restrict sparts = c->stars.parts;
165
  const struct engine *e = r->e;
166
167
  const struct unit_system *us = e->internal_units;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
168
  const struct cosmology *cosmo = e->cosmology;
169
  const struct feedback_props *feedback_props = e->feedback_props;
170
171
172
173
174
175
  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;
176
  int redo = 0, scount = 0;
177

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

181
182
  TIMER_TIC;

183
184
185
186
187
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID != e->nodeID)
    error("Running the star ghost on a foreign node!");
#endif

188
  /* Anything to do here? */
189
  if (c->stars.count == 0) return;
Loic Hausammann's avatar
Loic Hausammann committed
190
  if (!cell_is_active_stars(c, e)) return;
191
192
193

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

Matthieu Schaller's avatar
Matthieu Schaller committed
198
199
        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->stars.h_max);
200
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
201
    }
202
203
204
205
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
206
207
208
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
209
    if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
Loic Hausammann's avatar
Loic Hausammann committed
210
      error("Can't allocate memory for sid.");
211
    if ((h_0 = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
212
      error("Can't allocate memory for h_0.");
213
    if ((left = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
214
      error("Can't allocate memory for left.");
215
    if ((right = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
216
      error("Can't allocate memory for right.");
217
    for (int k = 0; k < c->stars.count; k++)
218
219
      if (spart_is_active(&sparts[k], e) &&
          feedback_is_active(&sparts[k], e->time, cosmo, with_cosmology)) {
220
        sid[scount] = k;
221
222
223
        h_0[scount] = sparts[k].h;
        left[scount] = 0.f;
        right[scount] = stars_h_max;
224
        ++scount;
225
226
227
      }

    /* While there are particles that need to be updated... */
228
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
229
230
231
232
233
234
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
235
      for (int i = 0; i < scount; i++) {
236
237
238
239
240
241

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
242
243
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
244
245
246
#endif

        /* Get some useful values */
247
        const float h_init = h_0[i];
248
249
250
        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);
251

252
253
254
        float h_new;
        int has_no_neighbours = 0;

255
        if (sp->density.wcount == 0.f) { /* No neighbours case */
256
257
258
259
260
261

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

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

263
264
265
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
266
          stars_end_density(sp, cosmo);
267
268

          /* Compute one step of the Newton-Raphson scheme */
269
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
270
          const float n_target = stars_eta_dim;
271
272
          const float f = n_sum - n_target;
          const float f_prime =
273
274
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
275

276
          /* Improve the bisection bounds */
277
278
279
280
          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);
281
282
283
284
285
286
287

#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

288
          /* Skip if h is already h_max and we don't have enough neighbours */
289
290
291
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
292

293
            stars_reset_feedback(sp);
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339

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

              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;
              if (with_cosmology) {
                star_age_end_of_step =
                    cosmology_get_delta_time_from_scale_factors(
                        cosmo, (double)sp->birth_scale_factor, cosmo->a);
              } else {
                star_age_end_of_step = (float)e->time - sp->birth_time;
              }

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

                /* Reset the feedback fields of the star particle */
                feedback_reset_feedback(sp, feedback_props);
              }
            } else {

              feedback_reset_feedback(sp, feedback_props);
            }
340

341
342
343
344
345
            /* Ok, we are done with this particle */
            continue;
          }

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

347
348
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
349
350
351
352
353
354
355
356
357
358
359
360

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

361
362
363
          /* 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);
364
365
366
367

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
368
369
370
371
372
373
        }

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

          /* Ok, correct then */
374
375
376
377
378
379
380
381
382
383
384
385
386
387

          /* 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;
          }
388
389

          /* If below the absolute maximum, try again */
390
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
391
392
393

            /* Flag for another round of fun */
            sid[redo] = sid[i];
394
395
396
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
397
398
399
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
400
            stars_init_spart(sp);
401
            feedback_init_spart(sp);
402
403
404
405

            /* Off we go ! */
            continue;

406
407
408
409
410
411
          } 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) {
412
413

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
414
            sp->h = stars_h_max;
415
416
417

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
418
              stars_spart_has_no_neighbours(sp, cosmo);
419
            }
420
421
422
423
424

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
425
426
427
          }
        }

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

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

433
        stars_reset_feedback(sp);
434

435
        /* Only do feedback if stars have a reasonable birth time */
436
        if (feedback_do_feedback(sp)) {
437

438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
          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;
453
          if (with_cosmology) {
454
            star_age_end_of_step = cosmology_get_delta_time_from_scale_factors(
455
456
                cosmo, sp->birth_scale_factor, (float)cosmo->a);
          } else {
457
            star_age_end_of_step = (float)e->time - sp->birth_time;
458
459
          }

460
461
462
463
464
465
466
467
468
469
          /* 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);
470
471
472
473
          } else {

            /* Reset the feedback fields of the star particle */
            feedback_reset_feedback(sp, feedback_props);
474
          }
475
476
477
478
        } else {

          /* Reset the feedback fields of the star particle */
          feedback_reset_feedback(sp, feedback_props);
479
        }
480
481
482
483
484
485
      }

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

      /* Re-set the counter for the next loop (potentially). */
486
487
      scount = redo;
      if (scount > 0) {
488
489
490
491
492

        /* 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
493
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
494
495
496
497
498
499
500
501

#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)
502
503
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
504
505
506
507
508
509

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

              /* Left or right? */
              if (l->t->ci == finger)
510
511
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
512
              else
513
514
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
515
516
517
518
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
519
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
520
                                                NULL, 1);
521
522
523
524
525
526

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

              /* Left or right? */
              if (l->t->ci == finger)
527
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
528
                                                  scount, l->t->cj, 1);
529
              else
530
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
531
                                                  scount, l->t->ci, 1);
532
533
534
535
536
537
            }
          }
        }
      }
    }

538
539
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
540
541
542
    }

    /* Be clean */
543
544
    free(left);
    free(right);
545
    free(sid);
546
    free(h_0);
547
548
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
549
550
  /* Update h_max */
  c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
551

552
  /* The ghost may not always be at the top level.
553
   * Therefore we need to update h_max between the super- and top-levels */
554
  if (c->stars.ghost) {
555
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
556
      atomic_max_d(&tmp->stars.h_max, h_max);
557
558
559
    }
  }

560
  if (timer) TIMER_TOC(timer_do_stars_ghost);
561
562
}

563
564
565
566
567
568
569
570
/**
 * @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 ?
 */
571
572
void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
                                         int timer) {
573
574
575
576
577
578

  struct bpart *restrict bparts = c->black_holes.parts;
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
  const float black_holes_h_max = e->hydro_properties->h_max;
  const float black_holes_h_min = e->hydro_properties->h_min;
579
  const float eps = e->black_holes_properties->h_tolerance;
580
  const float black_holes_eta_dim =
581
      pow_dimension(e->black_holes_properties->eta_neighbours);
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
  const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
  int redo = 0, bcount = 0;

  /* Running value of the maximal smoothing length */
  double h_max = c->black_holes.h_max;

  TIMER_TIC;

  /* Anything to do here? */
  if (c->black_holes.count == 0) return;
  if (!cell_is_active_black_holes(c, e)) return;

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) {
598
        runner_do_black_holes_density_ghost(r, c->progeny[k], 0);
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783

        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->black_holes.h_max);
      }
    }
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
    if ((sid = (int *)malloc(sizeof(int) * c->black_holes.count)) == NULL)
      error("Can't allocate memory for sid.");
    if ((h_0 = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
      error("Can't allocate memory for h_0.");
    if ((left = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
      error("Can't allocate memory for left.");
    if ((right = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
      error("Can't allocate memory for right.");
    for (int k = 0; k < c->black_holes.count; k++)
      if (bpart_is_active(&bparts[k], e)) {
        sid[bcount] = k;
        h_0[bcount] = bparts[k].h;
        left[bcount] = 0.f;
        right[bcount] = black_holes_h_max;
        ++bcount;
      }

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

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

      /* Loop over the remaining active parts in this cell. */
      for (int i = 0; i < bcount; i++) {

        /* Get a direct pointer on the part. */
        struct bpart *bp = &bparts[sid[i]];

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
        if (!bpart_is_active(bp, e))
          error("Ghost applied to inactive particle");
#endif

        /* Get some useful values */
        const float h_init = h_0[i];
        const float h_old = bp->h;
        const float h_old_dim = pow_dimension(h_old);
        const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);

        float h_new;
        int has_no_neighbours = 0;

        if (bp->density.wcount == 0.f) { /* No neighbours case */

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

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

        } else {

          /* Finish the density calculation */
          black_holes_end_density(bp, cosmo);

          /* Compute one step of the Newton-Raphson scheme */
          const float n_sum = bp->density.wcount * h_old_dim;
          const float n_target = black_holes_eta_dim;
          const float f = n_sum - n_target;
          const float f_prime =
              bp->density.wcount_dh * h_old_dim +
              hydro_dimension * bp->density.wcount * h_old_dim_minus_one;

          /* Improve the bisection bounds */
          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);

#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

          /* Skip if h is already h_max and we don't have enough neighbours */
          /* Same if we are below h_min */
          if (((bp->h >= black_holes_h_max) && (f < 0.f)) ||
              ((bp->h <= black_holes_h_min) && (f > 0.f))) {

            black_holes_reset_feedback(bp);

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

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

          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);

          /* 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, bp->id, h_init, h_old, h_new, f, f_prime, n_sum,
                n_target, left[i], right[i]);
          }

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

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

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

          /* Ok, correct then */

          /* 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 */
            bp->h = pow_inv_dimension(
                0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));

          } else {

            /* Normal case */
            bp->h = h_new;
          }

          /* If below the absolute maximum, try again */
          if (bp->h < black_holes_h_max && bp->h > black_holes_h_min) {

            /* Flag for another round of fun */
            sid[redo] = sid[i];
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
            redo += 1;

            /* Re-initialise everything */
            black_holes_init_bpart(bp);

            /* Off we go ! */
            continue;

          } else if (bp->h <= black_holes_h_min) {

            /* Ok, this particle is a lost cause... */
            bp->h = black_holes_h_min;

          } else if (bp->h >= black_holes_h_max) {

            /* Ok, this particle is a lost cause... */
            bp->h = black_holes_h_max;

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
              black_holes_bpart_has_no_neighbours(bp, cosmo);
            }

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

        /* We now have a particle whose smoothing length has converged */

784
785
        black_holes_reset_feedback(bp);

786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
        /* Check if h_max has increased */
        h_max = max(h_max, bp->h);
      }

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

      /* Re-set the counter for the next loop (potentially). */
      bcount = redo;
      if (bcount > 0) {

        /* Climb up the cell hierarchy. */
        for (struct cell *finger = c; finger != NULL; finger = finger->parent) {

          /* Run through this cell's density interactions. */
          for (struct link *l = finger->black_holes.density; l != NULL;
               l = l->next) {

#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)
              runner_doself_subset_branch_bh_density(r, finger, bparts, sid,
                                                     bcount);

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

              /* Left or right? */
              if (l->t->ci == finger)
                runner_dopair_subset_branch_bh_density(r, finger, bparts, sid,
                                                       bcount, l->t->cj);
              else
                runner_dopair_subset_branch_bh_density(r, finger, bparts, sid,
                                                       bcount, l->t->ci);
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
              runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
                                             NULL, 1);

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

              /* Left or right? */
              if (l->t->ci == finger)
                runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
                                               l->t->cj, 1);
              else
                runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
                                               l->t->ci, 1);
            }
          }
        }
      }
    }

    if (bcount) {
      error("Smoothing length failed to converge on %i particles.", bcount);
    }

    /* Be clean */
    free(left);
    free(right);
    free(sid);
    free(h_0);
  }

  /* Update h_max */
  c->black_holes.h_max = h_max;

  /* The ghost may not always be at the top level.
   * Therefore we need to update h_max between the super- and top-levels */
863
  if (c->black_holes.density_ghost) {
864
865
866
867
868
869
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
      atomic_max_d(&tmp->black_holes.h_max, h_max);
    }
  }

  if (timer) TIMER_TOC(timer_do_black_holes_ghost);
870
871
}

872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
/**
 * @brief Intermediate task after the BHs have done their swallowing step.
 * This is used to update the BH quantities if necessary.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
void runner_do_black_holes_swallow_ghost(struct runner *r, struct cell *c,
                                         int timer) {

  struct bpart *restrict bparts = c->black_holes.parts;
  const int count = c->black_holes.count;
  const struct engine *e = r->e;
  const int with_cosmology = e->policy & engine_policy_cosmology;

  TIMER_TIC;

  /* Anything to do here? */
  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_black_holes_swallow_ghost(r, c->progeny[k], 0);
  } else {

    /* Loop over the parts in this cell. */
    for (int i = 0; i < count; i++) {

      /* Get a direct pointer on the part. */
      struct bpart *bp = &bparts[i];

      if (bpart_is_active(bp, e)) {

        /* Get particle time-step */
        double dt;
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(bp->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(e->ti_current - 1, bp->time_bin);

          dt = cosmology_get_delta_time(e->cosmology, ti_begin,
                                        ti_begin + ti_step);
        } else {
          dt = get_timestep(bp->time_bin, e->time_base);
        }

        /* Compute variables required for the feedback loop */
        black_holes_prepare_feedback(bp, e->black_holes_properties,
                                     e->physical_constants, e->cosmology, dt);
      }
    }
  }

  if (timer) TIMER_TOC(timer_do_black_holes_ghost);
}

Tom Theuns's avatar
Tom Theuns committed
931
932
933
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
934
935
936
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
937
 */
938
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
939

940
941
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
942
943
944
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
945
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
946

947
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
948

949
  /* Anything to do here? */
950
  if (!cell_is_active_gravity(c, e)) return;
951

Tom Theuns's avatar
Tom Theuns committed
952
953
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
954
    for (int k = 0; k < 8; k++)
955
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
956
  } else {
957

958
959
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
960

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

964
      /* Is this part within the time step? */
965
      if (gpart_is_active(gp, e)) {
966
967
        external_gravity_acceleration(time, potential, constants, gp);
      }
968
    }
969
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
970

971
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
972
973
}

974
975
976
977
978
979
980
981
982
/**
 * @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) {

983
984
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
  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
1009
/**
1010
1011
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
1012
1013
1014
1015
1016
1017
1018
 *
 * @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) {

1019
  const struct engine *e = r->e;
1020
1021
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
1022
1023
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
1024
  const struct unit_system *us = e->internal_units;
1025
  const struct hydro_props *hydro_props = e->hydro_properties;
1026
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
1027
  const double time_base = e->time_base;
1028
  const integertime_t ti_current = e->ti_current;
1029
1030
1031
  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
1032
1033
1034

  TIMER_TIC;

1035
  /* Anything to do here? */
1036
  if (!cell_is_active_hydro(c, e)) return;
1037

Stefan Arridge's avatar
Stefan Arridge committed
1038
1039
1040
1041
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
1042
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
1043

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

1047
1048
1049
      /* 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
1050

1051
      if (part_is_active(p, e)) {
1052

1053
        double dt_cool, dt_therm;
1054
1055
1056
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
1057
1058
              get_integer_time_begin(ti_current - 1, p->time_bin);

1059
1060
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
1061
1062
1063
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

1064
1065
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
1066
          dt_therm = get_timestep(p->time_bin, time_base);
1067
        }
1068

1069
        /* Let's cool ! */
1070
1071
1072
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
1073
      }
Stefan Arridge's avatar
Stefan Arridge committed
1074
1075
1076
1077
1078
1079
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
1080
1081
1082
1083
1084
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

1085
  struct engine *e = r->e;
1086
  const struct cosmology *cosmo = e->cosmology;
1087
1088
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
1089
1090
1091
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
1092
  const int with_cosmology = (e->policy & engine_policy_cosmology);
1093
  const int with_feedback = (e->policy & engine_policy_feedback);
1094
1095
1096
  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;
1097
  const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
1098
1099
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
1100
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
1101
1102
1103

  TIMER_TIC;

1104
1105
1106
1107
1108
#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
1109
  /* Anything to do here? */
1110
  if (c->hydro.count == 0 || !cell_is_active_hydro(c, e)) {
1111
    star_formation_logger_log_inactive_cell(&c->stars.sfh);
1112
1113
    return;
  }
1114
1115

  /* Reset the SFR */
1116
  star_formation_logger_init(&c->stars.sfh);
Matthieu Schaller's avatar
Matthieu Schaller committed
1117
1118
1119
1120

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
1121
1122
1123
1124
1125
1126
1127
1128
      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 */
1129
        star_formation_logger_add(&c->stars.sfh, &cp->stars.sfh);
1130
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
1131
  } else {
1132
1133
1134
1135
1136
1137
1138
1139

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

1140
      /* Only work on active particles */
1141
1142
      if (part_is_active(p, e)) {

1143
1144
        /* 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
1145
1146
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
1147

1148
1149
1150
1151
1152
1153
          /* 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);
1154

1155
1156
1157
1158
1159
1160
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

1162
1163
1164
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
1165

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

1169
1170
1171
1172
1173
1174
1175
          /* 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);

1176
1177
1178
            /* Did we get a star? (Or did we run out of spare ones?) */
            if (sp != NULL) {

1179
1180
              /* message("We formed a star id=%lld cellID=%d", sp->id,
               * c->cellID); */
1181

1182
1183
              /* Copy the properties of the gas particle to the star particle */
              star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
Folkert Nobels's avatar
Folkert Nobels committed
1184
1185
                                             with_cosmology, phys_const,
                                             hydro_props, us, cooling);
1186
1187
1188

              /* Update the Star formation history */
              star_formation_logger_log_new_spart(sp, &c->stars.sfh);
1189
            }
1190
1191
1192
1193
1194
1195
1196
1197
1198
          }

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

Folkert Nobels's avatar
Folkert Nobels committed
1200
      } else { /* is active? */
1201

1202
        /* Check if the particle is not inhibited */
Matthieu Schaller's avatar