runner.c 141 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
#undef FUNCTION
123

124
125
/* 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
Tom Theuns's avatar
Tom Theuns committed
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
  TIMER_TIC;

174
175
176
177
178
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID != e->nodeID)
    error("Running the star ghost on a foreign node!");
#endif

179
  /* Anything to do here? */
180
  if (c->stars.count == 0) return;
Loic Hausammann's avatar
Loic Hausammann committed
181
  if (!cell_is_active_stars(c, e)) return;
182
183
184

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

Matthieu Schaller's avatar
Matthieu Schaller committed
189
190
        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->stars.h_max);
191
      }
Matthieu Schaller's avatar
Matthieu Schaller committed
192
    }
193
194
195
196
  } else {

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

    /* While there are particles that need to be updated... */
219
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
220
221
222
223
224
225
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
226
      for (int i = 0; i < scount; i++) {
227
228
229
230
231
232

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
233
234
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
235
236
237
#endif

        /* Get some useful values */
238
        const float h_init = h_0[i];
239
240
241
        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);
242

243
244
245
        float h_new;
        int has_no_neighbours = 0;

246
        if (sp->density.wcount == 0.f) { /* No neighbours case */
247
248
249
250
251
252

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

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

254
255
256
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
257
          stars_end_density(sp, cosmo);
258
259

          /* Compute one step of the Newton-Raphson scheme */
260
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
261
          const float n_target = stars_eta_dim;
262
263
          const float f = n_sum - n_target;
          const float f_prime =
264
265
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
266

267
          /* Improve the bisection bounds */
268
269
270
271
          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);
272
273
274
275
276
277
278

#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

279
          /* Skip if h is already h_max and we don't have enough neighbours */
280
281
282
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
283

284
            stars_reset_feedback(sp);
285
            feedback_reset_feedback(sp, feedback_props);
286

287
288
289
290
291
            /* Ok, we are done with this particle */
            continue;
          }

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

293
294
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
295
296
297
298
299
300
301
302
303
304
305
306

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

307
308
309
          /* 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);
310
311
312
313

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

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

          /* Ok, correct then */
320
321
322
323
324
325
326
327
328
329
330
331
332
333

          /* 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;
          }
334
335

          /* If below the absolute maximum, try again */
336
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
337
338
339

            /* Flag for another round of fun */
            sid[redo] = sid[i];
340
341
342
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
343
344
345
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
346
            stars_init_spart(sp);
347
            feedback_init_spart(sp);
348
349
350
351

            /* Off we go ! */
            continue;

352
353
354
355
356
357
          } 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) {
358
359

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
360
            sp->h = stars_h_max;
361
362
363

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
364
              stars_spart_has_no_neighbours(sp, cosmo);
365
            }
366
367
368
369
370

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

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

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

379
        stars_reset_feedback(sp);
380

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

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

406
407
408
409
410
411
412
413
414
415
          /* 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);
416
417
418
419
          } else {

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

          /* Reset the feedback fields of the star particle */
          feedback_reset_feedback(sp, feedback_props);
425
        }
426
427
428
429
430
431
      }

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

      /* Re-set the counter for the next loop (potentially). */
432
433
      scount = redo;
      if (scount > 0) {
434
435
436
437
438

        /* 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
439
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
440
441
442
443
444
445
446
447

#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)
448
449
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
450
451
452
453
454
455

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

              /* Left or right? */
              if (l->t->ci == finger)
456
457
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
458
              else
459
460
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
461
462
463
464
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
465
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
466
                                                NULL, 1);
467
468
469
470
471
472

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

              /* Left or right? */
              if (l->t->ci == finger)
473
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
474
                                                  scount, l->t->cj, 1);
475
              else
476
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
477
                                                  scount, l->t->ci, 1);
478
479
480
481
482
483
            }
          }
        }
      }
    }

484
485
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
486
487
488
    }

    /* Be clean */
489
490
    free(left);
    free(right);
491
    free(sid);
492
    free(h_0);
493
494
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
495
496
  /* Update h_max */
  c->stars.h_max = h_max;
Loic Hausammann's avatar
Loic Hausammann committed
497

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

506
  if (timer) TIMER_TOC(timer_do_stars_ghost);
507
508
}

509
510
511
512
513
514
515
516
517
518
519
520
/**
 * @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;
521
  const int with_cosmology = e->policy & engine_policy_cosmology;
522
523
524
  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;
525
  const float eps = e->black_holes_properties->h_tolerance;
526
  const float black_holes_eta_dim =
527
      pow_dimension(e->black_holes_properties->eta_neighbours);
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
635
636
637
638
639
  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))) {

640
641
642
643
644
645
646
647
648
649
650
651
652
            /* 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);
            }

653
            /* Compute variables required for the feedback loop */
654
655
656
            black_holes_prepare_feedback(bp, e->black_holes_properties,
                                         e->physical_constants, e->cosmology,
                                         dt);
657
658

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

752
753
754
755
756
757
758
759
760
761
762
763
764
        /* 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);
        }

765
        /* Compute variables required for the feedback loop */
766
767
        black_holes_prepare_feedback(bp, e->black_holes_properties,
                                     e->physical_constants, e->cosmology, dt);
768
769

        /* Reset quantities computed by the feedback loop */
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
850
851
852
        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);
853
854
}

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

864
865
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
866
867
868
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
869
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
870

871
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
872

873
  /* Anything to do here? */
874
  if (!cell_is_active_gravity(c, e)) return;
875

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

882
883
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
884

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

888
      /* Is this part within the time step? */
889
      if (gpart_is_active(gp, e)) {
890
891
        external_gravity_acceleration(time, potential, constants, gp);
      }
892
    }
893
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
894

895
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
896
897
}

898
899
900
901
902
903
904
905
906
/**
 * @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) {

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

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

  TIMER_TIC;

959
  /* Anything to do here? */
960
  if (!cell_is_active_hydro(c, e)) return;
961

Stefan Arridge's avatar
Stefan Arridge committed
962
963
964
965
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
966
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
967

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

971
972
973
      /* 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
974

975
      if (part_is_active(p, e)) {
976

977
        double dt_cool, dt_therm;
978
979
980
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
981
982
              get_integer_time_begin(ti_current - 1, p->time_bin);

983
984
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
985
986
987
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

988
989
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
990
          dt_therm = get_timestep(p->time_bin, time_base);
991
        }
992

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

For faster browsing, not all history is shown. View entire blame