runner.c 101 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"
Loikki's avatar
Loikki committed
56
#include "logger.h"
57
#include "minmax.h"
James Willis's avatar
James Willis committed
58
#include "runner_doiact_vec.h"
59
#include "scheduler.h"
60
#include "sort_part.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"
67
#include "timestep_limiter.h"
68
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
69

70
71
72
73
#define TASK_LOOP_DENSITY 0
#define TASK_LOOP_GRADIENT 1
#define TASK_LOOP_FORCE 2
#define TASK_LOOP_LIMITER 3
74

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

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

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

98
99
100
101
102
103
104
/* 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

105
/* Import the gravity loop functions. */
106
#include "runner_doiact_grav.h"
107

108
109
/* Import the stars density loop functions. */
#define FUNCTION density
Loic Hausammann's avatar
Loic Hausammann committed
110
#include "runner_doiact_stars.h"
111
112
113
114
#undef FUNCTION

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

/**
 * @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
126
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
127

128
  struct spart *restrict sparts = c->stars.parts;
129
130
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
Loic Hausammann's avatar
Loic Hausammann committed
131
  const struct stars_props *stars_properties = e->stars_properties;
Loic Hausammann's avatar
Loic Hausammann committed
132
  const float stars_h_max = stars_properties->h_max;
133
  const float eps = stars_properties->h_tolerance;
134
  const float stars_eta_dim = pow_dimension(stars_properties->eta_neighbours);
135
136
  const int max_smoothing_iter = stars_properties->max_smoothing_iterations;
  int redo = 0, scount = 0;
137
138
139
140

  TIMER_TIC;

  /* Anything to do here? */
Loic Hausammann's avatar
Loic Hausammann committed
141
  if (!cell_is_active_stars(c, e)) return;
142
143
144
145

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
Loic Hausammann's avatar
Loic Hausammann committed
146
      if (c->progeny[k] != NULL) runner_do_stars_ghost(r, c->progeny[k], 0);
147
148
149
150
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
151
    if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
Loic Hausammann's avatar
Loic Hausammann committed
152
      error("Can't allocate memory for sid.");
153
    for (int k = 0; k < c->stars.count; k++)
154
      if (spart_is_active(&sparts[k], e)) {
155
156
        sid[scount] = k;
        ++scount;
157
158
159
      }

    /* While there are particles that need to be updated... */
160
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
161
162
163
164
165
166
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
167
      for (int i = 0; i < scount; i++) {
168
169
170
171
172
173

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
174
175
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
176
177
178
179
180
181
182
183
184
#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;

185
        if (sp->density.wcount == 0.f) { /* No neighbours case */
186
187
188
189
190
191
192
193
194

          /* 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
195
          stars_end_density(sp, cosmo);
196
197

          /* Compute one step of the Newton-Raphson scheme */
198
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
199
          const float n_target = stars_eta_dim;
200
201
          const float f = n_sum - n_target;
          const float f_prime =
202
203
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224

          /* 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
225
          if (sp->h < stars_h_max) {
226
227
228
229
230
231

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

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
232
            stars_init_spart(sp);
233
234
235
236
237
238
239

            /* Off we go ! */
            continue;

          } else {

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
240
            sp->h = stars_h_max;
241
242
243

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
244
              stars_spart_has_no_neighbours(sp, cosmo);
245
246
247
248
            }
          }
        }

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

Loic Hausammann's avatar
Loic Hausammann committed
251
        /* Compute the stellar evolution  */
252
        stars_evolve_spart(sp, stars_properties, cosmo);
253
254
255
256
257
258
      }

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

      /* Re-set the counter for the next loop (potentially). */
259
260
      scount = redo;
      if (scount > 0) {
261
262
263
264
265

        /* 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
266
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
267
268
269
270
271
272
273
274

#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)
275
276
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
277
278
279
280
281
282

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

              /* Left or right? */
              if (l->t->ci == finger)
283
284
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
285
              else
286
287
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
288
289
290
291
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
292
293
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
294
295
296
297
298
299

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

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

311
312
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
313
314
315
316
317
318
    }

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

Matthieu Schaller's avatar
Matthieu Schaller committed
319
  if (timer) TIMER_TOC(timer_dostars_ghost);
320
321
}

Tom Theuns's avatar
Tom Theuns committed
322
323
324
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
325
326
327
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
328
 */
329
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
330

331
332
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
333
334
335
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
336
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
337

338
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
339

340
  /* Anything to do here? */
341
  if (!cell_is_active_gravity(c, e)) return;
342

Tom Theuns's avatar
Tom Theuns committed
343
344
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
345
    for (int k = 0; k < 8; k++)
346
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
347
  } else {
348

349
350
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
351

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

355
      /* Is this part within the time step? */
356
      if (gpart_is_active(gp, e)) {
357
358
        external_gravity_acceleration(time, potential, constants, gp);
      }
359
    }
360
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
361

362
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
363
364
}

365
366
367
368
369
370
371
372
373
/**
 * @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) {

374
375
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
  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
400
/**
401
402
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
403
404
405
406
407
408
409
 *
 * @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) {

410
  const struct engine *e = r->e;
411
412
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
413
414
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
415
  const struct unit_system *us = e->internal_units;
416
  const struct hydro_props *hydro_props = e->hydro_properties;
417
  const double time_base = e->time_base;
418
  const integertime_t ti_current = e->ti_current;
419
420
421
  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
422
423
424

  TIMER_TIC;

425
  /* Anything to do here? */
426
  if (!cell_is_active_hydro(c, e)) return;
427

Stefan Arridge's avatar
Stefan Arridge committed
428
429
430
431
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
432
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
433

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

437
438
439
      /* 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
440

441
      if (part_is_active(p, e)) {
442

443
        double dt_cool, dt_therm;
444
445
446
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
447
448
              get_integer_time_begin(ti_current - 1, p->time_bin);

449
450
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
451
452
453
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

454
455
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
456
          dt_therm = get_timestep(p->time_bin, time_base);
457
        }
458

459
        /* Let's cool ! */
460
461
        cooling_cool_part(constants, us, cosmo, hydro_props, cooling_func, p,
                          xp, dt_cool, dt_therm);
462
      }
Stefan Arridge's avatar
Stefan Arridge committed
463
464
465
466
467
468
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
469
470
471
472
473
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

474
  struct engine *e = r->e;
475
  const struct cosmology *cosmo = e->cosmology;
476
477
478
  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
479
480
481
482
483
484
485
486
487
488
489

  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 {
490
491
492
493
494
495
496
497
498
499

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

500
501
        const float rho = hydro_get_physical_density(p, cosmo);

502
        // MATTHIEU: Temporary star-formation law
503
        // Do not use this at home.
504
        if (rho > 1.7e7 && e->step > 2) {
505
          message("Removing particle id=%lld rho=%e", p->id, rho);
506

507
          struct spart *sp = cell_convert_part_to_spart(e, c, p, xp);
508

509
          /* Did we run out of fresh particles? */
510
          if (sp == NULL) continue;
511

512
513
          /* Set everything to a valid state */
          stars_init_spart(sp);
514
515
516
        }
      }
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
517
518
519
520
521
  }

  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
522
523
524
525
526
527
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
528
void runner_do_sort_ascending(struct entry *sort, int N) {
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

  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
583
        }
584
585
586
587
588
589
590
591
592
593
594
595
      } 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
596
    }
597
598
599
  }
}

600
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
601
602
603
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
604
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
605
 */
606
607
608
609
610
611
612
613
614
#define RUNNER_CHECK_SORTS(TYPE)                                               \
  void runner_check_sorts_##TYPE(struct cell *c, int flags) {                  \
                                                                               \
    if (flags & ~c->TYPE.sorted) error("Inconsistent sort flags (downward)!"); \
    if (c->split)                                                              \
      for (int k = 0; k < 8; k++)                                              \
        if (c->progeny[k] != NULL && c->progeny[k]->TYPE.count > 0)            \
          runner_check_sorts_##TYPE(c->progeny[k], c->TYPE.sorted);            \
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
615
#else
Loic Hausammann's avatar
Loic Hausammann committed
616
617
618
619
#define RUNNER_CHECK_SORTS(TYPE)                                       \
  void runner_check_sorts_##TYPE(struct cell *c, int flags) {          \
    error("Calling debugging code without debugging flag activated."); \
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
620
#endif
Loic Hausammann's avatar
Loic Hausammann committed
621
622
623

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
624

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

  struct entry *fingers[8];
640
641
642
  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
643
  float buff[8];
644

645
646
  TIMER_TIC;

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

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

660
661
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
662
  runner_check_sorts_hydro(c, c->hydro.sorted);
663
664

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

671
672
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
673
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
674
#endif
675

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

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

    /* Fill in the gaps within the progeny. */
689
    float dx_max_sort = 0.0f;
690
    float dx_max_sort_old = 0.0f;
691
    for (int k = 0; k < 8; k++) {
692
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
693
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
694
        runner_do_hydro_sort(r, c->progeny[k], flags,
695
696
697
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
698
699
700
        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);
701
      }
702
    }
703
704
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
705
706

    /* Loop over the 13 different sort arrays. */
707
    for (int j = 0; j < 13; j++) {
708
709
710
711
712

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

      /* Init the particle index offsets. */
713
      int off[8];
714
715
      off[0] = 0;
      for (int k = 1; k < 8; k++)
716
        if (c->progeny[k - 1] != NULL)
717
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
718
719
720
721
        else
          off[k] = off[k - 1];

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

      /* Sort the buffer. */
734
735
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
736
          if (buff[inds[k]] < buff[inds[i]]) {
737
            int temp_i = inds[i];
738
739
740
741
742
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
743
      struct entry *finger = c->hydro.sort[j];
744
      for (int ind = 0; ind < count; ind++) {
745
746
747
748
749
750
751
752
753
754

        /* 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. */
755
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
756
          int temp_i = inds[k - 1];
757
758
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
759
        }
760

761
762
763
      } /* Merge. */

      /* Add a sentinel. */
764
765
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
766
767

      /* Mark as sorted. */
768
      atomic_or(&c->hydro.sorted, 1 << j);
769
770
771
772
773
774
775
776

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

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

      /* 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;
        }
791
      }
792
793
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
794
795
    }

796
    /* Fill the sort array. */
797
    for (int k = 0; k < count; k++) {
798
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
799
      for (int j = 0; j < 13; j++)
800
        if (flags & (1 << j)) {
801
802
803
804
          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];
805
        }
806
    }
807
808

    /* Add the sentinel and sort. */
809
    for (int j = 0; j < 13; j++)
810
      if (flags & (1 << j)) {
811
812
813
814
        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);
815
816
817
      }
  }

818
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
819
  /* Verify the sorting. */
820
  for (int j = 0; j < 13; j++) {
821
    if (!(flags & (1 << j))) continue;
822
    struct entry *finger = c->hydro.sort[j];
823
    for (int k = 1; k < count; k++) {
824
825
826
827
828
      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
829

830
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
831
  runner_check_sorts_hydro(c, flags);
832
833

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

841
  /* Clear the cell's sort flags. */
842
843
844
  c->hydro.do_sort = 0;
  c->hydro.do_sub_sort = 0;
  c->hydro.requires_sorts = 0;
845

846
847
848
  if (clock) TIMER_TOC(timer_dosort);
}

Loic Hausammann's avatar
Loic Hausammann committed
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
/**
 * @brief Sort the stars particles in the given cell along all cardinal
 * directions.
 *
 * @param r The #runner.
 * @param c The #cell.
 * @param flags Cell flag.
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
 */
void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {

  struct entry *fingers[8];
  const int count = c->stars.count;
  struct spart *sparts = c->stars.parts;
  float buff[8];

  TIMER_TIC;

  /* We need to do the local sorts plus whatever was requested further up. */
  flags |= c->stars.do_sort;
  if (cleanup) {
    c->stars.sorted = 0;
  } else {
    flags &= ~c->stars.sorted;
  }
  if (flags == 0 && !c->stars.do_sub_sort) return;

  /* Check that the particles have been moved to the current time */
  if (flags && !cell_are_spart_drifted(c, r->e))
    error("Sorting un-drifted cell c->nodeID=%d", c->nodeID);

#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
  runner_check_sorts_stars(c, c->stars.sorted);

  /* Make sure the sort flags are consistent (upward). */
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
    if (finger->stars.sorted & ~c->stars.sorted)
      error("Inconsistent sort flags (upward).");
  }

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

  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
    if ((flags & (1 << j)) && c->stars.sort[j] == NULL) {
      if ((c->stars.sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                                     (count + 1))) == NULL)
        error("Failed to allocate sort memory.");
    }
  }

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

    /* Fill in the gaps within the progeny. */
    float dx_max_sort = 0.0f;
    float dx_max_sort_old = 0.0f;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL && c->progeny[k]->stars.count > 0) {
        /* Only propagate cleanup if the progeny is stale. */
        runner_do_stars_sort(r, c->progeny[k], flags,
Loic Hausammann's avatar
Loic Hausammann committed
919
                             cleanup && (c->progeny[k]->stars.dx_max_sort_old >
Loic Hausammann's avatar
Loic Hausammann committed
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
        dx_max_sort = max(dx_max_sort, c->progeny[k]->stars.dx_max_sort);
        dx_max_sort_old =
            max(dx_max_sort_old, c->progeny[k]->stars.dx_max_sort_old);
      }
    }
    c->stars.dx_max_sort = dx_max_sort;
    c->stars.dx_max_sort_old = dx_max_sort_old;

    /* Loop over the 13 different sort arrays. */
    for (int j = 0; j < 13; j++) {

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

      /* Init the particle index offsets. */
      int off[8];
      off[0] = 0;
      for (int k = 1; k < 8; k++)
        if (c->progeny[k - 1] != NULL)
          off[k] = off[k - 1] + c->progeny[k - 1]->stars.count;
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
      int inds[8];
      for (int k = 0; k < 8; k++) {
        inds[k] = k;
        if (c->progeny[k] != NULL && c->progeny[k]->stars.count > 0) {
          fingers[k] = c->progeny[k]->stars.sort[j];
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
          if (buff[inds[k]] < buff[inds[i]]) {
            int temp_i = inds[i];
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
      struct entry *finger = c->stars.sort[j];
      for (int ind = 0; ind < count; ind++) {

        /* 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. */
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
          int temp_i = inds[k - 1];
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
        }

      } /* Merge. */

      /* Add a sentinel. */
      c->stars.sort[j][count].d = FLT_MAX;
      c->stars.sort[j][count].i = 0;

      /* Mark as sorted. */
      atomic_or(&c->stars.sorted, 1 << j);

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

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