runner.c 138 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
79
80
81
#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
82
#define TASK_LOOP_FEEDBACK 4
83

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

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

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

107
108
109
110
111
112
113
/* 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

114
/* Import the gravity loop functions. */
115
#include "runner_doiact_grav.h"
116

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

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

131
132
133
134
135
136
137
138
139
140
141
142
143
144
/* 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

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

145
146
147
148
149
150
151
152
/**
 * @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
153
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
154

155
  struct spart *restrict sparts = c->stars.parts;
156
  const struct engine *e = r->e;
157
158
  const struct unit_system *us = e->internal_units;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
159
  const struct cosmology *cosmo = e->cosmology;
160
  const struct feedback_props *feedback_props = e->feedback_props;
161
162
163
164
165
166
  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;
167
  int redo = 0, scount = 0;
168

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

172
173
174
  TIMER_TIC;

  /* Anything to do here? */
175
  if (c->stars.count == 0) return;
Loic Hausammann's avatar
Loic Hausammann committed
176
  if (!cell_is_active_stars(c, e)) return;
177
178
179

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

Matthieu Schaller's avatar
Matthieu Schaller committed
184
185
        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->stars.h_max);
186
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
187
    }
188
189
190
191
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
192
193
194
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
195
    if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
Loic Hausammann's avatar
Loic Hausammann committed
196
      error("Can't allocate memory for sid.");
197
    if ((h_0 = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
198
      error("Can't allocate memory for h_0.");
199
    if ((left = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
200
      error("Can't allocate memory for left.");
201
    if ((right = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
202
      error("Can't allocate memory for right.");
203
    for (int k = 0; k < c->stars.count; k++)
204
205
      if (spart_is_active(&sparts[k], e) &&
          feedback_is_active(&sparts[k], e->time, cosmo, with_cosmology)) {
206
        sid[scount] = k;
207
208
209
        h_0[scount] = sparts[k].h;
        left[scount] = 0.f;
        right[scount] = stars_h_max;
210
        ++scount;
211
212
213
      }

    /* While there are particles that need to be updated... */
214
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
215
216
217
218
219
220
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
221
      for (int i = 0; i < scount; i++) {
222
223
224
225
226
227

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
228
229
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
230
231
232
#endif

        /* Get some useful values */
233
        const float h_init = h_0[i];
234
235
236
        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);
237

238
239
240
        float h_new;
        int has_no_neighbours = 0;

241
        if (sp->density.wcount == 0.f) { /* No neighbours case */
242
243
244
245
246
247

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

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

249
250
251
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
252
          stars_end_density(sp, cosmo);
253
254

          /* Compute one step of the Newton-Raphson scheme */
255
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
256
          const float n_target = stars_eta_dim;
257
258
          const float f = n_sum - n_target;
          const float f_prime =
259
260
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
261

262
          /* Improve the bisection bounds */
263
264
265
266
          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);
267
268
269
270
271
272
273

#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

274
          /* Skip if h is already h_max and we don't have enough neighbours */
275
276
277
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
278

279
            stars_reset_feedback(sp);
280
            feedback_reset_feedback(sp, feedback_props);
281

282
283
284
285
286
            /* Ok, we are done with this particle */
            continue;
          }

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

288
289
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
290
291
292
293
294
295
296
297
298
299
300
301

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

302
303
304
          /* 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);
305
306
307
308

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

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

          /* Ok, correct then */
315
316
317
318
319
320
321
322
323
324
325
326
327
328

          /* 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;
          }
329
330

          /* If below the absolute maximum, try again */
331
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
332
333
334

            /* Flag for another round of fun */
            sid[redo] = sid[i];
335
336
337
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
338
339
340
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
341
            stars_init_spart(sp);
342
            feedback_init_spart(sp);
343
344
345
346

            /* Off we go ! */
            continue;

347
348
349
350
351
352
          } 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) {
353
354

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
355
            sp->h = stars_h_max;
356
357
358

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
359
              stars_spart_has_no_neighbours(sp, cosmo);
360
            }
361
362
363
364
365

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

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

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

374
        stars_reset_feedback(sp);
375

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

379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
          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;
394
          if (with_cosmology) {
395
            star_age_end_of_step = cosmology_get_delta_time_from_scale_factors(
396
397
                cosmo, sp->birth_scale_factor, (float)cosmo->a);
          } else {
398
            star_age_end_of_step = (float)e->time - sp->birth_time;
399
400
          }

401
402
403
404
405
406
407
408
409
410
          /* 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);
411
412
413
414
          } else {

            /* Reset the feedback fields of the star particle */
            feedback_reset_feedback(sp, feedback_props);
415
          }
416
417
418
419
        } else {

          /* Reset the feedback fields of the star particle */
          feedback_reset_feedback(sp, feedback_props);
420
        }
421
422
423
424
425
426
      }

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

      /* Re-set the counter for the next loop (potentially). */
427
428
      scount = redo;
      if (scount > 0) {
429
430
431
432
433

        /* 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
434
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
435
436
437
438
439
440
441
442

#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)
443
444
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
445
446
447
448
449
450

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

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

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
460
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
461
                                                NULL, 1);
462
463
464
465
466
467

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

              /* Left or right? */
              if (l->t->ci == finger)
468
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
469
                                                  scount, l->t->cj, 1);
470
              else
471
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
472
                                                  scount, l->t->ci, 1);
473
474
475
476
477
478
            }
          }
        }
      }
    }

479
480
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
481
482
483
    }

    /* Be clean */
484
485
    free(left);
    free(right);
486
    free(sid);
487
    free(h_0);
488
489
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
490
491
  /* Update h_max */
  c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
492

493
  /* The ghost may not always be at the top level.
494
   * Therefore we need to update h_max between the super- and top-levels */
495
  if (c->stars.ghost) {
496
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
497
      atomic_max_d(&tmp->stars.h_max, h_max);
498
499
500
    }
  }

501
  if (timer) TIMER_TOC(timer_do_stars_ghost);
502
503
}

504
505
506
507
508
509
510
511
512
513
514
515
/**
 * @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 ?
 */
void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {

  struct bpart *restrict bparts = c->black_holes.parts;
  const struct engine *e = r->e;
516
  const int with_cosmology = e->policy & engine_policy_cosmology;
517
518
519
  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;
520
  const float eps = e->black_holes_properties->h_tolerance;
521
  const float black_holes_eta_dim =
522
      pow_dimension(e->black_holes_properties->eta_neighbours);
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
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
  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) {
        runner_do_black_holes_ghost(r, c->progeny[k], 0);

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

635
636
637
638
639
640
641
642
643
644
645
646
647
            /* 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);
            }

648
            /* Compute variables required for the feedback loop */
649
650
651
            black_holes_prepare_feedback(bp, e->black_holes_properties,
                                         e->physical_constants, e->cosmology,
                                         dt);
652
653

            /* Reset quantities computed by the feedback loop */
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
            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 */

        /* Check if h_max has increased */
        h_max = max(h_max, bp->h);

747
748
749
750
751
752
753
754
755
756
757
758
759
        /* 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);
        }

760
        /* Compute variables required for the feedback loop */
761
762
        black_holes_prepare_feedback(bp, e->black_holes_properties,
                                     e->physical_constants, e->cosmology, dt);
763
764

        /* Reset quantities computed by the feedback loop */
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
        black_holes_reset_feedback(bp);
      }

      /* 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 */
  if (c->black_holes.ghost) {
    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);
}

Tom Theuns's avatar
Tom Theuns committed
850
851
852
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
853
854
855
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
856
 */
857
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
858

859
860
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
861
862
863
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
864
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
865

866
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
867

868
  /* Anything to do here? */
869
  if (!cell_is_active_gravity(c, e)) return;
870

Tom Theuns's avatar
Tom Theuns committed
871
872
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
873
    for (int k = 0; k < 8; k++)
874
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
875
  } else {
876

877
878
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
879

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

883
      /* Is this part within the time step? */
884
      if (gpart_is_active(gp, e)) {
885
886
        external_gravity_acceleration(time, potential, constants, gp);
      }
887
    }
888
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
889

890
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
891
892
}

893
894
895
896
897
898
899
900
901
/**
 * @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) {

902
903
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
  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
928
/**
929
930
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
931
932
933
934
935
936
937
 *
 * @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) {

938
  const struct engine *e = r->e;
939
940
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
941
942
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
943
  const struct unit_system *us = e->internal_units;
944
  const struct hydro_props *hydro_props = e->hydro_properties;
945
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
946
  const double time_base = e->time_base;
947
  const integertime_t ti_current = e->ti_current;
948
949
950
  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
951
952
953

  TIMER_TIC;

954
  /* Anything to do here? */
955
  if (!cell_is_active_hydro(c, e)) return;
956

Stefan Arridge's avatar
Stefan Arridge committed
957
958
959
960
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
961
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
962

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

966
967
968
      /* 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
969

970
      if (part_is_active(p, e)) {
971

972
        double dt_cool, dt_therm;
973
974
975
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
976
977
              get_integer_time_begin(ti_current - 1, p->time_bin);

978
979
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
980
981
982
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

983
984
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
985
          dt_therm = get_timestep(p->time_bin, time_base);
986
        }
987

988
        /* Let's cool ! */
989
990
991
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
992
      }
Stefan Arridge's avatar
Stefan Arridge committed
993
994
995
996
997
998
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
999
1000
/**
 *
For faster browsing, not all history is shown. View entire blame