runner.c 94.4 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"
Folkert Nobels's avatar
Folkert Nobels committed
63
#include "starformation.h"
64
#include "stars.h"
65
66
#include "task.h"
#include "timers.h"
67
#include "timestep.h"
Folkert Nobels's avatar
Folkert Nobels committed
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
/* Import the gravity loop functions. */
99
#include "runner_doiact_grav.h"
100

101
102
/* Import the stars density loop functions. */
#define FUNCTION density
Loic Hausammann's avatar
Loic Hausammann committed
103
#include "runner_doiact_stars.h"
104
105
106
107
#undef FUNCTION

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

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

121
  struct spart *restrict sparts = c->stars.parts;
122
123
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
Loic Hausammann's avatar
Loic Hausammann committed
124
  const struct stars_props *stars_properties = e->stars_properties;
Loic Hausammann's avatar
Loic Hausammann committed
125
  const float stars_h_max = stars_properties->h_max;
126
  const float eps = stars_properties->h_tolerance;
127
  const float stars_eta_dim = pow_dimension(stars_properties->eta_neighbours);
128
129
  const int max_smoothing_iter = stars_properties->max_smoothing_iterations;
  int redo = 0, scount = 0;
130
131
132
133

  TIMER_TIC;

  /* Anything to do here? */
Loic Hausammann's avatar
Loic Hausammann committed
134
  if (!cell_is_active_stars(c, e)) return;
135
136
137
138

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
Loic Hausammann's avatar
Loic Hausammann committed
139
      if (c->progeny[k] != NULL) runner_do_stars_ghost(r, c->progeny[k], 0);
140
141
142
143
  } else {

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

    /* While there are particles that need to be updated... */
153
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
154
155
156
157
158
159
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
160
      for (int i = 0; i < scount; i++) {
161
162
163
164
165
166

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
167
168
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
169
170
171
172
173
174
175
176
177
#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;

178
        if (sp->density.wcount == 0.f) { /* No neighbours case */
179
180
181
182
183
184
185
186
187

          /* 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
188
          stars_end_density(sp, cosmo);
189
190

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

          /* 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
218
          if (sp->h < stars_h_max) {
219
220
221
222
223
224

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

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
225
            stars_init_spart(sp);
226
227
228
229
230
231
232

            /* Off we go ! */
            continue;

          } else {

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
233
            sp->h = stars_h_max;
234
235
236

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
237
              stars_spart_has_no_neighbours(sp, cosmo);
238
239
240
241
            }
          }
        }

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

Loic Hausammann's avatar
Loic Hausammann committed
244
        /* Compute the stellar evolution  */
245
        stars_evolve_spart(sp, stars_properties, cosmo);
246
247
248
249
250
251
      }

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

      /* Re-set the counter for the next loop (potentially). */
252
253
      scount = redo;
      if (scount > 0) {
254
255
256
257
258

        /* 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
259
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
260
261
262
263
264
265
266
267

#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)
268
269
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
270
271
272
273
274
275

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

              /* Left or right? */
              if (l->t->ci == finger)
276
277
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
278
              else
279
280
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
281
282
283
284
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
285
286
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
287
288
289
290
291
292

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

              /* Left or right? */
              if (l->t->ci == finger)
293
294
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
295
              else
296
297
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
298
299
300
301
302
303
            }
          }
        }
      }
    }

304
305
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
306
307
308
309
310
311
    }

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

Matthieu Schaller's avatar
Matthieu Schaller committed
312
  if (timer) TIMER_TOC(timer_dostars_ghost);
313
314
}

Tom Theuns's avatar
Tom Theuns committed
315
316
317
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
318
319
320
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
321
 */
322
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
323

324
325
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
326
327
328
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
329
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
330

331
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
332

333
  /* Anything to do here? */
334
  if (!cell_is_active_gravity(c, e)) return;
335

Tom Theuns's avatar
Tom Theuns committed
336
337
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
338
    for (int k = 0; k < 8; k++)
339
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
340
  } else {
341

342
343
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
344

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

348
      /* Is this part within the time step? */
349
      if (gpart_is_active(gp, e)) {
350
351
        external_gravity_acceleration(time, potential, constants, gp);
      }
352
    }
353
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
354

355
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
356
357
}

358
359
360
361
362
363
364
365
366
/**
 * @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) {

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

403
  const struct engine *e = r->e;
404
405
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
406
407
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
408
  const struct unit_system *us = e->internal_units;
409
  const struct hydro_props *hydro_props = e->hydro_properties;
410
  const double time_base = e->time_base;
411
  const integertime_t ti_current = e->ti_current;
412
413
414
  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
415
416
417

  TIMER_TIC;

418
  /* Anything to do here? */
419
  if (!cell_is_active_hydro(c, e)) return;
420

Stefan Arridge's avatar
Stefan Arridge committed
421
422
423
424
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
425
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
426

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

430
431
432
      /* 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
433

434
      if (part_is_active(p, e)) {
435

436
        double dt_cool, dt_therm;
437
438
439
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
440
441
              get_integer_time_begin(ti_current - 1, p->time_bin);

442
443
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
444
445
446
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

447
448
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
449
          dt_therm = get_timestep(p->time_bin, time_base);
450
        }
451

452
        /* Let's cool ! */
453
454
        cooling_cool_part(constants, us, cosmo, hydro_props, cooling_func, p,
                          xp, dt_cool, dt_therm);
455
      }
Stefan Arridge's avatar
Stefan Arridge committed
456
457
458
459
460
461
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
462
463
464
465
466
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

467
  struct engine *e = r->e;
468
  const struct cosmology *cosmo = e->cosmology;
469
  const struct star_formation *starform = e->star_formation;
Folkert Nobels's avatar
Folkert Nobels committed
470
  const struct phys_const *constants = e->physical_constants;
471
472
473
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
474
  const int with_cosmology = (e->policy & engine_policy_cosmology);
475
476
477
  const struct hydro_props *restrict hydro_props = e->hydro_properties;
  const struct unit_system *restrict us = e->internal_units;
  struct cooling_function_data *restrict cooling = e->cooling_func;
478
479
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
480
481
482
483
484
485
486
487
488
489
490

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

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

501
502
        /* Calculate the time step of the current particle for cosmo and no
         * cosmo*/
503
504
505
506
507
508
509
510
511
512
513
514
515
        double dt_star;
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(ti_current - 1, p->time_bin);

          dt_star =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

Folkert Nobels's avatar
Folkert Nobels committed
516
517
        // const float rho = hydro_get_physical_density(p, cosmo);
        if (star_formation_convert_to_star(e, starform, p, xp, constants, cosmo,
Matthieu Schaller's avatar
Matthieu Schaller committed
518
519
                                           hydro_props, us, cooling, dt_star,
                                           with_cosmology)) {
520
          /* Convert your particle to a star */
521
          struct spart *sp = cell_convert_part_to_spart(e, c, p, xp);
522

523
          /* Copy the properties of the gas particle to the star particle */
524
          star_formation_copy_properties(e, p, xp, sp, starform, constants,
Folkert Nobels's avatar
Folkert Nobels committed
525
                                         cosmo, with_cosmology);
Folkert Nobels's avatar
Folkert Nobels committed
526
        }
527
528
      }
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
529
530
531
532
533
  }

  if (timer) TIMER_TOC(timer_do_star_formation);
}

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

  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
595
        }
596
597
598
599
600
601
602
603
604
605
606
607
      } 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
608
    }
609
610
611
  }
}

612
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
613
614
615
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
616
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
617
 */
618
619
620
621
622
623
624
625
626
#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
627
#else
Loic Hausammann's avatar
Loic Hausammann committed
628
629
630
631
#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
632
#endif
Loic Hausammann's avatar
Loic Hausammann committed
633
634
635

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
636

Pedro Gonnet's avatar
Pedro Gonnet committed
637
638
639
640
641
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
642
 * @param flags Cell flag.
643
644
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
645
646
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
647
 */
Loic Hausammann's avatar
Loic Hausammann committed
648
649
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
650
651

  struct entry *fingers[8];
652
653
654
  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
655
  float buff[8];
656

657
658
  TIMER_TIC;

659
  /* We need to do the local sorts plus whatever was requested further up. */
660
  flags |= c->hydro.do_sort;
661
  if (cleanup) {
662
    c->hydro.sorted = 0;
663
  } else {
664
    flags &= ~c->hydro.sorted;
665
  }
666
  if (flags == 0 && !c->hydro.do_sub_sort) return;
667
668

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

672
673
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
674
  runner_check_sorts_hydro(c, c->hydro.sorted);
675
676

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
677
678
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
679
680
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
681
  }
682

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

688
689
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
690
691
692
    if ((flags & (1 << j)) && c->hydro.sort[j] == NULL) {
      if ((c->hydro.sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                                     (count + 1))) == NULL)
693
694
        error("Failed to allocate sort memory.");
    }
695
696
697
698
699
700
  }

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

    /* Fill in the gaps within the progeny. */
701
    float dx_max_sort = 0.0f;
702
    float dx_max_sort_old = 0.0f;
703
    for (int k = 0; k < 8; k++) {
704
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
705
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
706
        runner_do_hydro_sort(r, c->progeny[k], flags,
707
708
709
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
710
711
712
        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);
713
      }
714
    }
715
716
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
717
718

    /* Loop over the 13 different sort arrays. */
719
    for (int j = 0; j < 13; j++) {
720
721
722
723
724

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

      /* Init the particle index offsets. */
725
      int off[8];
726
727
      off[0] = 0;
      for (int k = 1; k < 8; k++)
728
        if (c->progeny[k - 1] != NULL)
729
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
730
731
732
733
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
734
      int inds[8];
735
      for (int k = 0; k < 8; k++) {
736
        inds[k] = k;
737
738
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
739
740
741
742
743
744
745
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
746
747
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
748
          if (buff[inds[k]] < buff[inds[i]]) {
749
            int temp_i = inds[i];
750
751
752
753
754
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
755
      struct entry *finger = c->hydro.sort[j];
756
      for (int ind = 0; ind < count; ind++) {
757
758
759
760
761
762
763
764
765
766

        /* 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. */
767
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
768
          int temp_i = inds[k - 1];
769
770
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
771
        }
772

773
774
775
      } /* Merge. */

      /* Add a sentinel. */
776
777
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
778
779

      /* Mark as sorted. */
780
      atomic_or(&c->hydro.sorted, 1 << j);
781
782
783
784
785
786
787
788

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

789
    /* Reset the sort distance */
790
    if (c->hydro.sorted == 0) {
791
792
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
793
        error("Have non-NULL xparts in foreign cell");
794
#endif
795
796
797
798
799
800
801
802

      /* 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;
        }
803
      }
804
805
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
806
807
    }

808
    /* Fill the sort array. */
809
    for (int k = 0; k < count; k++) {
810
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
811
      for (int j = 0; j < 13; j++)
812
        if (flags & (1 << j)) {
813
814
815
816
          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];
817
        }
818
    }
819
820

    /* Add the sentinel and sort. */
821
    for (int j = 0; j < 13; j++)
822
      if (flags & (1 << j)) {
823
824
825
826
        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);
827
828
829
      }
  }

830
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
831
  /* Verify the sorting. */
832
  for (int j = 0; j < 13; j++) {
833
    if (!(flags & (1 << j))) continue;
834
    struct entry *finger = c->hydro.sort[j];
835
    for (int k = 1; k < count; k++) {
836
837
838
839
840
      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
841

842
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
843
  runner_check_sorts_hydro(c, flags);
844
845

  /* Make sure the sort flags are consistent (upward). */
Pedro Gonnet's avatar
Pedro Gonnet committed
846
847
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
848
849
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags.");
850
  }
851
#endif
852

853
  /* Clear the cell's sort flags. */
854
855
856
  c->hydro.do_sort = 0;
  c->hydro.do_sub_sort = 0;
  c->hydro.requires_sorts = 0;
857

858
859
860
  if (clock) TIMER_TOC(timer_dosort);
}

Loic Hausammann's avatar
Loic Hausammann committed
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
919
920
921
922
923
924
925
926
927
928
929
930
/**
 * @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
931
                             cleanup && (c->progeny[k]->stars.dx_max_sort_old >
Loic Hausammann's avatar
Loic Hausammann committed
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;
For faster browsing, not all history is shown. View entire blame