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

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

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

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

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

Pedro Gonnet's avatar
Pedro Gonnet committed
40
/* Local headers. */
41
#include "active.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
42
#include "approx_math.h"
43
#include "atomic.h"
44
#include "cell.h"
45
#include "chemistry.h"
46
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
47
#include "cooling.h"
48
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
49
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
50
#include "engine.h"
51
#include "error.h"
52
53
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
54
#include "hydro_properties.h"
55
#include "kick.h"
56
#include "minmax.h"
James Willis's avatar
James Willis committed
57
#include "runner_doiact_vec.h"
58
#include "scheduler.h"
59
#include "sort_part.h"
60
#include "sourceterms.h"
61
#include "space.h"
62
#include "space_getsid.h"
63
#include "stars.h"
64
65
#include "task.h"
#include "timers.h"
66
#include "timestep.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
67

68
69
70
71
#define TASK_LOOP_DENSITY 0
#define TASK_LOOP_GRADIENT 1
#define TASK_LOOP_FORCE 2
#define TASK_LOOP_LIMITER 3
72

73
/* Import the density loop functions. */
74
#define FUNCTION density
75
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
76
#include "runner_doiact.h"
77
78
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
79

80
/* Import the gradient loop functions (if required). */
81
82
#ifdef EXTRA_HYDRO_LOOP
#define FUNCTION gradient
83
#define FUNCTION_TASK_LOOP TASK_LOOP_GRADIENT
84
#include "runner_doiact.h"
85
86
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
87
88
#endif

89
/* Import the force loop functions. */
90
#define FUNCTION force
91
#define FUNCTION_TASK_LOOP TASK_LOOP_FORCE
92
#include "runner_doiact.h"
93
94
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
95

96
/* Import the gravity loop functions. */
97
#include "runner_doiact_grav.h"
98

Loic Hausammann's avatar
Loic Hausammann committed
99
100
/* Import the stars loop functions. */
#include "runner_doiact_stars.h"
101

Tom Theuns's avatar
Tom Theuns committed
102
/**
Tom Theuns's avatar
Tom Theuns committed
103
 * @brief Perform source terms
Tom Theuns's avatar
Tom Theuns committed
104
105
106
107
108
109
 *
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
 */
void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
110
  const int count = c->hydro.count;
111
  const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
Tom Theuns's avatar
Tom Theuns committed
112
  const double cell_width[3] = {c->width[0], c->width[1], c->width[2]};
Tom Theuns's avatar
Tom Theuns committed
113
  struct sourceterms *sourceterms = r->e->sourceterms;
114
  const int dimen = 3;
Tom Theuns's avatar
Tom Theuns committed
115
116
117
118
119
120
121

  TIMER_TIC;

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_sourceterms(r, c->progeny[k], 0);
122
  } else {
Tom Theuns's avatar
Tom Theuns committed
123

124
    if (count > 0) {
Tom Theuns's avatar
Tom Theuns committed
125

126
127
128
129
130
131
      /* do sourceterms in this cell? */
      const int incell =
          sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
      if (incell == 1) {
        sourceterms_apply(r, sourceterms, c);
      }
Tom Theuns's avatar
Tom Theuns committed
132
133
    }
  }
Tom Theuns's avatar
Tom Theuns committed
134
135
136
137

  if (timer) TIMER_TOC(timer_dosource);
}

138
139
140
141
142
143
144
145
/**
 * @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
146
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
147

148
  struct spart *restrict sparts = c->stars.parts;
149
150
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
Loic Hausammann's avatar
Loic Hausammann committed
151
  const struct stars_props *stars_properties = e->stars_properties;
Loic Hausammann's avatar
Loic Hausammann committed
152
  const float stars_h_max = stars_properties->h_max;
153
  const float eps = stars_properties->h_tolerance;
154
  const float stars_eta_dim = pow_dimension(stars_properties->eta_neighbours);
155
156
  const int max_smoothing_iter = stars_properties->max_smoothing_iterations;
  int redo = 0, scount = 0;
157
158
159
160

  TIMER_TIC;

  /* Anything to do here? */
Loic Hausammann's avatar
Loic Hausammann committed
161
  if (!cell_is_active_stars(c, e)) return;
162
163
164
165

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
Loic Hausammann's avatar
Loic Hausammann committed
166
      if (c->progeny[k] != NULL) runner_do_stars_ghost(r, c->progeny[k], 0);
167
168
169
170
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
171
    if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
Loic Hausammann's avatar
Loic Hausammann committed
172
      error("Can't allocate memory for sid.");
173
    for (int k = 0; k < c->stars.count; k++)
174
      if (spart_is_active(&sparts[k], e)) {
175
176
        sid[scount] = k;
        ++scount;
177
178
179
      }

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

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

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

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

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

        /* Get some useful values */
        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);
        float h_new;
        int has_no_neighbours = 0;

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

          /* 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 */
Loic Hausammann's avatar
Loic Hausammann committed
215
          stars_end_density(sp, cosmo);
216
217

          /* Compute one step of the Newton-Raphson scheme */
218
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
219
          const float n_target = stars_eta_dim;
220
221
          const float f = n_sum - n_target;
          const float f_prime =
222
223
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244

          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
#ifdef SWIFT_DEBUG_CHECKS
          if ((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old))
            error(
                "Smoothing length correction not going in the right direction");
#endif

          /* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
          h_new = min(h_new, 2.f * h_old);
          h_new = max(h_new, 0.5f * h_old);
        }

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

          /* Ok, correct then */
          sp->h = h_new;

          /* If below the absolute maximum, try again */
Loic Hausammann's avatar
Loic Hausammann committed
245
          if (sp->h < stars_h_max) {
246
247
248
249
250
251

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

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
252
            stars_init_spart(sp);
253
254
255
256
257
258
259

            /* Off we go ! */
            continue;

          } else {

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
260
            sp->h = stars_h_max;
261
262
263

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
264
              stars_spart_has_no_neighbours(sp, cosmo);
265
266
267
268
            }
          }
        }

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

Loic Hausammann's avatar
Loic Hausammann committed
271
        /* Compute the stellar evolution  */
272
        stars_evolve_spart(sp, stars_properties, cosmo);
273
274
275
276
277
278
      }

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

      /* Re-set the counter for the next loop (potentially). */
279
280
      scount = redo;
      if (scount > 0) {
281
282
283
284
285

        /* 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
286
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
287
288
289
290
291
292
293
294

#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)
295
296
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
297
298
299
300
301
302

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

              /* Left or right? */
              if (l->t->ci == finger)
303
304
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
305
              else
306
307
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
308
309
310
311
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
312
313
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
314
315
316
317
318
319

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

              /* Left or right? */
              if (l->t->ci == finger)
320
321
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
322
              else
323
324
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
325
326
327
328
329
330
            }
          }
        }
      }
    }

331
332
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
333
334
335
336
337
338
    }

    /* Be clean */
    free(sid);
  }

Loic Hausammann's avatar
Loic Hausammann committed
339
  if (timer) TIMER_TOC(timer_do_stars_ghost);
340
341
}

Tom Theuns's avatar
Tom Theuns committed
342
343
344
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
345
346
347
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
348
 */
349
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
350

351
352
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
353
354
355
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
356
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
357

358
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
359

360
  /* Anything to do here? */
361
  if (!cell_is_active_gravity(c, e)) return;
362

Tom Theuns's avatar
Tom Theuns committed
363
364
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
365
    for (int k = 0; k < 8; k++)
366
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
367
  } else {
368

369
370
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
371

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

375
      /* Is this part within the time step? */
376
      if (gpart_is_active(gp, e)) {
377
378
        external_gravity_acceleration(time, potential, constants, gp);
      }
379
    }
380
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
381

382
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
383
384
}

385
386
387
388
389
390
391
392
393
/**
 * @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) {

394
395
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
  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
420
/**
421
422
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
423
424
425
426
427
428
429
 *
 * @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) {

430
  const struct engine *e = r->e;
431
432
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
433
434
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
435
  const struct unit_system *us = e->internal_units;
436
  const struct hydro_props *hydro_props = e->hydro_properties;
437
  const double time_base = e->time_base;
438
  const integertime_t ti_current = e->ti_current;
439
440
441
  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
442
443
444

  TIMER_TIC;

445
  /* Anything to do here? */
446
  if (!cell_is_active_hydro(c, e)) return;
447

Stefan Arridge's avatar
Stefan Arridge committed
448
449
450
451
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
452
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
453

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

457
458
459
      /* 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
460

461
      if (part_is_active(p, e)) {
462

463
        double dt_cool, dt_therm;
464
465
466
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
467
468
              get_integer_time_begin(ti_current - 1, p->time_bin);

469
470
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
471
472
473
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

474
475
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
476
          dt_therm = get_timestep(p->time_bin, time_base);
477
        }
478

479
        /* Let's cool ! */
480
481
        cooling_cool_part(constants, us, cosmo, hydro_props, cooling_func, p,
                          xp, dt_cool, dt_therm);
482
      }
Stefan Arridge's avatar
Stefan Arridge committed
483
484
485
486
487
488
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
489
490
491
492
493
494
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

  const struct engine *e = r->e;
495
496
497
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
498
499
500
501
502
503
504
505
506
507
508

  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_star_formation(r, c->progeny[k], 0);
  } else {
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525

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

      if (part_is_active(p, e)) {

        // MATTHIEU: Temporary star-formation law
        if (p->rho > 1.5e7 && e->step > 2) {
          message("Removing particle id=%lld rho=%e", p->id, p->rho);
          cell_convert_part_to_gpart(e, c, p, xp);
        }
      }
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
526
527
528
529
530
  }

  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
531
532
533
534
535
536
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
537
void runner_do_sort_ascending(struct entry *sort, int N) {
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

  struct {
    short int lo, hi;
  } qstack[10];
  int qpos, i, j, lo, hi, imin;
  struct entry temp;
  float pivot;

  /* Sort parts in cell_i in decreasing order with quicksort */
  qstack[0].lo = 0;
  qstack[0].hi = N - 1;
  qpos = 0;
  while (qpos >= 0) {
    lo = qstack[qpos].lo;
    hi = qstack[qpos].hi;
    qpos -= 1;
    if (hi - lo < 15) {
      for (i = lo; i < hi; i++) {
        imin = i;
        for (j = i + 1; j <= hi; j++)
          if (sort[j].d < sort[imin].d) imin = j;
        if (imin != i) {
          temp = sort[imin];
          sort[imin] = sort[i];
          sort[i] = temp;
        }
      }
    } else {
      pivot = sort[(lo + hi) / 2].d;
      i = lo;
      j = hi;
      while (i <= j) {
        while (sort[i].d < pivot) i++;
        while (sort[j].d > pivot) j--;
        if (i <= j) {
          if (i < j) {
            temp = sort[i];
            sort[i] = sort[j];
            sort[j] = temp;
          }
          i += 1;
          j -= 1;
        }
      }
      if (j > (lo + hi) / 2) {
        if (lo < j) {
          qpos += 1;
          qstack[qpos].lo = lo;
          qstack[qpos].hi = j;
        }
        if (i < hi) {
          qpos += 1;
          qstack[qpos].lo = i;
          qstack[qpos].hi = hi;
Pedro Gonnet's avatar
Pedro Gonnet committed
592
        }
593
594
595
596
597
598
599
600
601
602
603
604
      } else {
        if (i < hi) {
          qpos += 1;
          qstack[qpos].lo = i;
          qstack[qpos].hi = hi;
        }
        if (lo < j) {
          qpos += 1;
          qstack[qpos].lo = lo;
          qstack[qpos].hi = j;
        }
      }
Pedro Gonnet's avatar
Pedro Gonnet committed
605
    }
606
607
608
  }
}

Matthieu Schaller's avatar
Matthieu Schaller committed
609
610
611
612
613
614
615
616
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
 * Debugging function.
 *
 * @param c The #cell to check.
 * @param flags The sorting flags to check.
 */
617
void runner_check_sorts(struct cell *c, int flags) {
Matthieu Schaller's avatar
Matthieu Schaller committed
618
619

#ifdef SWIFT_DEBUG_CHECKS
620
  if (flags & ~c->hydro.sorted) error("Inconsistent sort flags (downward)!");
621
622
  if (c->split)
    for (int k = 0; k < 8; k++)
623
624
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0)
        runner_check_sorts(c->progeny[k], c->hydro.sorted);
Matthieu Schaller's avatar
Matthieu Schaller committed
625
626
627
#else
  error("Calling debugging code without debugging flag activated.");
#endif
628
629
}

Pedro Gonnet's avatar
Pedro Gonnet committed
630
631
632
633
634
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
635
 * @param flags Cell flag.
636
637
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
638
639
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
640
 */
641
642
void runner_do_sort(struct runner *r, struct cell *c, int flags, int cleanup,
                    int clock) {
643
644

  struct entry *fingers[8];
645
646
647
  const int count = c->hydro.count;
  const struct part *parts = c->hydro.parts;
  struct xpart *xparts = c->hydro.xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
648
  float buff[8];
649

650
651
  TIMER_TIC;

652
  /* We need to do the local sorts plus whatever was requested further up. */
653
  flags |= c->hydro.do_sort;
654
  if (cleanup) {
655
    c->hydro.sorted = 0;
656
  } else {
657
    flags &= ~c->hydro.sorted;
658
  }
659
  if (flags == 0 && !c->hydro.do_sub_sort) return;
660
661

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

665
666
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
667
  runner_check_sorts(c, c->hydro.sorted);
668
669

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
670
671
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
672
673
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
674
  }
675

676
677
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
678
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
679
#endif
680

681
682
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
683
684
685
    if ((flags & (1 << j)) && c->hydro.sort[j] == NULL) {
      if ((c->hydro.sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                                     (count + 1))) == NULL)
686
687
        error("Failed to allocate sort memory.");
    }
688
689
690
691
692
693
  }

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

    /* Fill in the gaps within the progeny. */
694
    float dx_max_sort = 0.0f;
695
    float dx_max_sort_old = 0.0f;
696
    for (int k = 0; k < 8; k++) {
697
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
698
699
        /* Only propagate cleanup if the progeny is stale. */
        runner_do_sort(r, c->progeny[k], flags,
700
                       cleanup && (c->progeny[k]->hydro.dx_max_sort >
701
702
                                   space_maxreldx * c->progeny[k]->dmin),
                       0);
703
704
705
        dx_max_sort = max(dx_max_sort, c->progeny[k]->hydro.dx_max_sort);
        dx_max_sort_old =
            max(dx_max_sort_old, c->progeny[k]->hydro.dx_max_sort_old);
706
      }
707
    }
708
709
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
710
711

    /* Loop over the 13 different sort arrays. */
712
    for (int j = 0; j < 13; j++) {
713
714
715
716
717

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

      /* Init the particle index offsets. */
718
      int off[8];
719
720
      off[0] = 0;
      for (int k = 1; k < 8; k++)
721
        if (c->progeny[k - 1] != NULL)
722
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
723
724
725
726
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
727
      int inds[8];
728
      for (int k = 0; k < 8; k++) {
729
        inds[k] = k;
730
731
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
732
733
734
735
736
737
738
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
739
740
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
741
          if (buff[inds[k]] < buff[inds[i]]) {
742
            int temp_i = inds[i];
743
744
745
746
747
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
748
      struct entry *finger = c->hydro.sort[j];
749
      for (int ind = 0; ind < count; ind++) {
750
751
752
753
754
755
756
757
758
759

        /* Copy the minimum into the new sort array. */
        finger[ind].d = buff[inds[0]];
        finger[ind].i = fingers[inds[0]]->i + off[inds[0]];

        /* Update the buffer. */
        fingers[inds[0]] += 1;
        buff[inds[0]] = fingers[inds[0]]->d;

        /* Find the smallest entry. */
760
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
761
          int temp_i = inds[k - 1];
762
763
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
764
        }
765

766
767
768
      } /* Merge. */

      /* Add a sentinel. */
769
770
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
771
772

      /* Mark as sorted. */
773
      atomic_or(&c->hydro.sorted, 1 << j);
774
775
776
777
778
779
780
781

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

782
    /* Reset the sort distance */
783
    if (c->hydro.sorted == 0) {
784
785
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
786
        error("Have non-NULL xparts in foreign cell");
787
#endif
788
789
790
791
792
793
794
795

      /* And the individual sort distances if we are a local cell */
      if (xparts != NULL) {
        for (int k = 0; k < count; k++) {
          xparts[k].x_diff_sort[0] = 0.0f;
          xparts[k].x_diff_sort[1] = 0.0f;
          xparts[k].x_diff_sort[2] = 0.0f;
        }
796
      }
797
798
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
799
800
    }

801
    /* Fill the sort array. */
802
    for (int k = 0; k < count; k++) {
803
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
804
      for (int j = 0; j < 13; j++)
805
        if (flags & (1 << j)) {
806
807
808
809
          c->hydro.sort[j][k].i = k;
          c->hydro.sort[j][k].d = px[0] * runner_shift[j][0] +
                                  px[1] * runner_shift[j][1] +
                                  px[2] * runner_shift[j][2];
810
        }
811
    }
812
813

    /* Add the sentinel and sort. */
814
    for (int j = 0; j < 13; j++)
815
      if (flags & (1 << j)) {
816
817
818
819
        c->hydro.sort[j][count].d = FLT_MAX;
        c->hydro.sort[j][count].i = 0;
        runner_do_sort_ascending(c->hydro.sort[j], count);
        atomic_or(&c->hydro.sorted, 1 << j);
820
821
822
      }
  }

823
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
824
  /* Verify the sorting. */
825
  for (int j = 0; j < 13; j++) {
826
    if (!(flags & (1 << j))) continue;
827
    struct entry *finger = c->hydro.sort[j];
828
    for (int k = 1; k < count; k++) {
829
830
831
832
833
      if (finger[k].d < finger[k - 1].d)
        error("Sorting failed, ascending array.");
      if (finger[k].i >= count) error("Sorting failed, indices borked.");
    }
  }
Pedro Gonnet's avatar
Pedro Gonnet committed
834

835
836
837
838
  /* Make sure the sort flags are consistent (downward). */
  runner_check_sorts(c, flags);

  /* Make sure the sort flags are consistent (upward). */
Pedro Gonnet's avatar
Pedro Gonnet committed
839
840
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
841
842
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags.");
843
  }
844
#endif
845

846
  /* Clear the cell's sort flags. */
847
848
849
  c->hydro.do_sort = 0;
  c->hydro.do_sub_sort = 0;
  c->hydro.requires_sorts = 0;
850

851
852
853
  if (clock) TIMER_TOC(timer_dosort);
}

854
/**
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
 * @brief Initialize the multipoles before the gravity calculation.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer 1 if the time is to be recorded.
 */
void runner_do_init_grav(struct runner *r, struct cell *c, int timer) {

  const struct engine *e = r->e;

  TIMER_TIC;

#ifdef SWIFT_DEBUG_CHECKS
  if (!(e->policy & engine_policy_self_gravity))
    error("Grav-init task called outside of self-gravity calculation");
#endif

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

  /* Reset the gravity acceleration tensors */
876
  gravity_field_tensors_init(&c->grav.multipole->pot, e->ti_current);
877
878
879
880
881
882
883
884
885
886
887

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) runner_do_init_grav(r, c->progeny[k], 0);
    }
  }

  if (timer) TIMER_TOC(timer_init_grav);
}

888
/**
889
890
891
892
893
 * @brief Intermediate task after the gradient loop that does final operations
 * on the gradient quantities and optionally slope limits the gradients
 *
 * @param r The runner thread.
 * @param c The cell.
894
 * @param timer Are we timing this ?
895
 */
896
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
897

898
#ifdef EXTRA_HYDRO_LOOP
899

900
901
902
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
  const int count = c->hydro.count;
903
  const struct engine *e = r->e;
904
  const integertime_t ti_end = e->ti_current;
Josh Borrow's avatar
Josh Borrow committed
905
  const int with_cosmology = (e->policy & engine_policy_cosmology);
906
  const double time_base = e->time_base;
907
  const struct cosmology *cosmo = e->cosmology;
Josh Borrow's avatar
Josh Borrow committed
908
  const struct hydro_props *hydro_props = e->hydro_properties;
909

910
911
  TIMER_TIC;

912
  /* Anything to do here? */
913
  if (!cell_is_active_hydro(c, e)) return;
914

915
916
917
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
918
      if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k], 0);
919
920
921
922
923
924
925
  } else {

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

      /* Get a direct pointer on the part. */
      struct part *restrict p = &parts[i];
926
      struct xpart *restrict xp = &xparts[i];
927

928
      if (part_is_active(p, e)) {
929

930
        /* Finish the gradient calculation */
931
        hydro_end_gradient(p);
932
933
934

        /* As of here, particle force variables will be set. */

Josh Borrow's avatar
Josh Borrow committed
935
936
937
        /* Calculate the time-step for passing to hydro_prepare_force.
         * This is the physical time between the start and end of the time-step
         * without any scale-factor powers. */
938
939
940
941
942
943
944
        double dt_alpha;
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          dt_alpha = cosmology_get_delta_time(cosmo, ti_end - ti_step, ti_end);
        } else {
          dt_alpha = get_timestep(p->time_bin, time_base);
        }
945

946
        /* Compute variables required for the force loop */
Josh Borrow's avatar
Josh Borrow committed
947
        hydro_prepare_force(p, xp, cosmo, hydro_props, dt_alpha);
948
949
950
951
952
953

        /* The particle force values are now set.  Do _NOT_
           try to read any particle density variables! */

        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
954
955
956
      }
    }
  }
957

958
959
  if (timer) TIMER_TOC(timer_do_extra_ghost);

960
961
#else
  error("SWIFT was not compiled with the extra hydro loop activated.");
962
#endif
963
}
964

965
/**
966
967
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
968
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
969
 * @param r The runner thread.
970
 * @param c The cell.
971
 * @param timer Are we timing this ?
972
 */
973
void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
974

975
976
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
977
  const struct engine *e = r->e;
Bert Vandenbroucke's avatar