runner.c 104 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 "entropy_floor.h"
52
#include "error.h"
53
54
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
55
#include "hydro_properties.h"
56
#include "kick.h"
Loikki's avatar
Loikki committed
57
#include "logger.h"
58
#include "minmax.h"
James Willis's avatar
James Willis committed
59
#include "runner_doiact_vec.h"
60
#include "scheduler.h"
61
#include "sort_part.h"
62
#include "space.h"
63
#include "space_getsid.h"
64
#include "star_formation.h"
65
#include "stars.h"
66
67
#include "task.h"
#include "timers.h"
68
#include "timestep.h"
69
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
70
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
71

72
73
74
75
#define TASK_LOOP_DENSITY 0
#define TASK_LOOP_GRADIENT 1
#define TASK_LOOP_FORCE 2
#define TASK_LOOP_LIMITER 3
76

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

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

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

100
101
102
103
104
105
106
/* 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

107
/* Import the gravity loop functions. */
108
#include "runner_doiact_grav.h"
109

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

/* Import the stars feedback loop functions. */
#define FUNCTION feedback
Loic Hausammann's avatar
Loic Hausammann committed
119
#include "runner_doiact_stars.h"
120
#undef FUNCTION
121
122
123
124
125
126
127
128
129

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

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

  TIMER_TIC;

  /* Anything to do here? */
Loic Hausammann's avatar
Loic Hausammann committed
145
  if (!cell_is_active_stars(c, e)) return;
146
147
148
149

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

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

    /* While there are particles that need to be updated... */
164
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
165
166
167
168
169
170
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
171
      for (int i = 0; i < scount; i++) {
172
173
174
175
176
177

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

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

189
        if (sp->density.wcount == 0.f) { /* No neighbours case */
190
191
192
193
194
195
196
197
198

          /* 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
199
          stars_end_density(sp, cosmo);
200
201

          /* Compute one step of the Newton-Raphson scheme */
202
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
203
          const float n_target = stars_eta_dim;
204
205
          const float f = n_sum - n_target;
          const float f_prime =
206
207
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
208

209
210
211
          /* Skip if h is already h_max and we don't have enough neighbours */
          if ((sp->h >= stars_h_max) && (f < 0.f)) {

212
            stars_reset_acceleration(sp);
213
214
215
216
217
            /* Ok, we are done with this particle */
            continue;
          }

          /* Normal case: Use Newton-Raphson to get a better value of h */
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
          /* 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
238
          if (sp->h < stars_h_max) {
239
240
241
242
243
244

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

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
245
            stars_init_spart(sp);
246
247
248
249
250
251
252

            /* Off we go ! */
            continue;

          } else {

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
253
            sp->h = stars_h_max;
254
255
256

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
257
              stars_spart_has_no_neighbours(sp, cosmo);
258
259
260
261
            }
          }
        }

262
        /* We now have a particle whose smoothing length has converged */
263
        stars_reset_acceleration(sp);
264

Loic Hausammann's avatar
Loic Hausammann committed
265
        /* Compute the stellar evolution  */
266
        stars_evolve_spart(sp, stars_properties, cosmo);
267
268
269
270
271
272
      }

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

      /* Re-set the counter for the next loop (potentially). */
273
274
      scount = redo;
      if (scount > 0) {
275
276
277
278
279

        /* 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
280
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
281
282
283
284
285
286
287
288

#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)
289
290
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
291
292
293
294
295
296

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

              /* Left or right? */
              if (l->t->ci == finger)
297
298
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
299
              else
300
301
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
302
303
304
305
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
306
307
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
308
309
310
311
312
313

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

              /* Left or right? */
              if (l->t->ci == finger)
314
315
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
316
              else
317
318
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
319
320
321
322
323
324
            }
          }
        }
      }
    }

325
326
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
327
328
329
330
331
332
    }

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

Matthieu Schaller's avatar
Matthieu Schaller committed
333
  if (timer) TIMER_TOC(timer_dostars_ghost);
334
335
}

Tom Theuns's avatar
Tom Theuns committed
336
337
338
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
339
340
341
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
342
 */
343
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
344

345
346
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
347
348
349
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
350
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
351

352
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
353

354
  /* Anything to do here? */
355
  if (!cell_is_active_gravity(c, e)) return;
356

Tom Theuns's avatar
Tom Theuns committed
357
358
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
359
    for (int k = 0; k < 8; k++)
360
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
361
  } else {
362

363
364
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
365

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

369
      /* Is this part within the time step? */
370
      if (gpart_is_active(gp, e)) {
371
372
        external_gravity_acceleration(time, potential, constants, gp);
      }
373
    }
374
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
375

376
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
377
378
}

379
380
381
382
383
384
385
386
387
/**
 * @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) {

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

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

  TIMER_TIC;

439
  /* Anything to do here? */
440
  if (!cell_is_active_hydro(c, e)) return;
441

Stefan Arridge's avatar
Stefan Arridge committed
442
443
444
445
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
446
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
447

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

451
452
453
      /* 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
454

455
      if (part_is_active(p, e)) {
456

457
        double dt_cool, dt_therm;
458
459
460
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
461
462
              get_integer_time_begin(ti_current - 1, p->time_bin);

463
464
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
465
466
467
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

468
469
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
470
          dt_therm = get_timestep(p->time_bin, time_base);
471
        }
472

473
        /* Let's cool ! */
474
475
        cooling_cool_part(constants, us, cosmo, hydro_props, cooling_func, p,
                          xp, dt_cool, dt_therm);
476
      }
Stefan Arridge's avatar
Stefan Arridge committed
477
478
479
480
481
482
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
483
484
485
486
487
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

488
  struct engine *e = r->e;
489
  const struct cosmology *cosmo = e->cosmology;
490
491
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
492
493
494
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
495
  const int with_cosmology = (e->policy & engine_policy_cosmology);
496
497
498
  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;
499
500
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
501
502
503
504
505
506
507
508
509
510
511

  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 {
512
513
514
515
516
517
518
519

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

520
      /* Only work on active particles */
521
522
      if (part_is_active(p, e)) {

523
524
525
        /* Is this particle star forming? */
        if (star_formation_is_star_forming(p, xp, sf_props, phys_const, cosmo,
                                           hydro_props, us, cooling)) {
526

527
528
529
530
531
532
          /* Time-step size for this particle */
          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);
533

534
535
536
537
538
539
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

541
542
543
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
544

545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
          /* Are we forming a star particle from this SF rate? */
          if (star_formation_should_convert_to_star(p, xp, sf_props, e,
                                                    dt_star)) {

            /* Convert the gas particle to a star particle */
            struct spart *sp = cell_convert_part_to_spart(e, c, p, xp);

            /* Copy the properties of the gas particle to the star particle */
            star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
                                           with_cosmology);
          }

        } else { /* Are we not star-forming? */

          /* Update the particle to flag it as not star-forming */
          star_formation_update_part_not_SFR(p, xp, e, sf_props,
                                             with_cosmology);

        } /* Not Star-forming? */
      }   /* is active? */
    }     /* Loop over particles */
Matthieu Schaller's avatar
Matthieu Schaller committed
566
567
568
569
570
  }

  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
571
572
573
574
575
576
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
577
void runner_do_sort_ascending(struct entry *sort, int N) {
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631

  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
632
        }
633
634
635
636
637
638
639
640
641
642
643
644
      } 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
645
    }
646
647
648
  }
}

649
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
650
651
652
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
653
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
654
 */
655
656
657
658
659
660
661
662
663
#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
664
#else
Loic Hausammann's avatar
Loic Hausammann committed
665
666
667
668
#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
669
#endif
Loic Hausammann's avatar
Loic Hausammann committed
670
671
672

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
673

Pedro Gonnet's avatar
Pedro Gonnet committed
674
675
676
677
678
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
679
 * @param flags Cell flag.
680
681
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
682
683
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
684
 */
Loic Hausammann's avatar
Loic Hausammann committed
685
686
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
687
688

  struct entry *fingers[8];
689
690
691
  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
692
  float buff[8];
693

694
695
  TIMER_TIC;

696
  /* We need to do the local sorts plus whatever was requested further up. */
697
  flags |= c->hydro.do_sort;
698
  if (cleanup) {
699
    c->hydro.sorted = 0;
700
  } else {
701
    flags &= ~c->hydro.sorted;
702
  }
703
  if (flags == 0 && !c->hydro.do_sub_sort) return;
704
705

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

709
710
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
711
  runner_check_sorts_hydro(c, c->hydro.sorted);
712
713

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
714
715
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
716
717
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
718
  }
719

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

725
726
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
727
728
729
    if ((flags & (1 << j)) && c->hydro.sort[j] == NULL) {
      if ((c->hydro.sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                                     (count + 1))) == NULL)
730
731
        error("Failed to allocate sort memory.");
    }
732
733
734
735
736
737
  }

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

    /* Fill in the gaps within the progeny. */
738
    float dx_max_sort = 0.0f;
739
    float dx_max_sort_old = 0.0f;
740
    for (int k = 0; k < 8; k++) {
741
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
742
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
743
        runner_do_hydro_sort(r, c->progeny[k], flags,
744
745
746
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
747
748
749
        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);
750
      }
751
    }
752
753
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
754
755

    /* Loop over the 13 different sort arrays. */
756
    for (int j = 0; j < 13; j++) {
757
758
759
760
761

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

      /* Init the particle index offsets. */
762
      int off[8];
763
764
      off[0] = 0;
      for (int k = 1; k < 8; k++)
765
        if (c->progeny[k - 1] != NULL)
766
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
767
768
769
770
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
771
      int inds[8];
772
      for (int k = 0; k < 8; k++) {
773
        inds[k] = k;
774
775
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
776
777
778
779
780
781
782
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
783
784
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
785
          if (buff[inds[k]] < buff[inds[i]]) {
786
            int temp_i = inds[i];
787
788
789
790
791
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
792
      struct entry *finger = c->hydro.sort[j];
793
      for (int ind = 0; ind < count; ind++) {
794
795
796
797
798
799
800
801
802
803

        /* 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. */
804
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
805
          int temp_i = inds[k - 1];
806
807
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
808
        }
809

810
811
812
      } /* Merge. */

      /* Add a sentinel. */
813
814
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
815
816

      /* Mark as sorted. */
817
      atomic_or(&c->hydro.sorted, 1 << j);
818
819
820
821
822
823
824
825

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

826
    /* Reset the sort distance */
827
    if (c->hydro.sorted == 0) {
828
829
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
830
        error("Have non-NULL xparts in foreign cell");
831
#endif
832
833
834
835
836
837
838
839

      /* 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;
        }
840
      }
841
842
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
843
844
    }

845
    /* Fill the sort array. */
846
    for (int k = 0; k < count; k++) {
847
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
848
      for (int j = 0; j < 13; j++)
849
        if (flags & (1 << j)) {
850
851
852
853
          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];
854
        }
855
    }
856
857

    /* Add the sentinel and sort. */
858
    for (int j = 0; j < 13; j++)
859
      if (flags & (1 << j)) {
860
861
862
863
        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);
864
865
866
      }
  }

867
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
868
  /* Verify the sorting. */
869
  for (int j = 0; j < 13; j++) {
870
    if (!(flags & (1 << j))) continue;
871
    struct entry *finger = c->hydro.sort[j];
872
    for (int k = 1; k < count; k++) {
873
874
875
876
877
      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
878

879
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
880
  runner_check_sorts_hydro(c, flags);
881
882

  /* Make sure the sort flags are consistent (upward). */
Pedro Gonnet's avatar
Pedro Gonnet committed
883
884
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
885
886
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags.");
887
  }
888
#endif
889

890
  /* Clear the cell's sort flags. */
891
892
893
  c->hydro.do_sort = 0;
  c->hydro.do_sub_sort = 0;
  c->hydro.requires_sorts = 0;
894

895
896
897
  if (clock) TIMER_TOC(timer_dosort);
}

Loic Hausammann's avatar
Loic Hausammann committed
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
/**
 * @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 */
930
  if (flags && !cell_are_spart_drifted(c, r->e)) {
Loic Hausammann's avatar
Loic Hausammann committed
931
    error("Sorting un-drifted cell c->nodeID=%d", c->nodeID);
932
  }
Loic Hausammann's avatar
Loic Hausammann committed
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

#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. */
968
969
970
971
        const int cleanup_prog =
            cleanup && (c->progeny[k]->stars.dx_max_sort_old >
                        space_maxreldx * c->progeny[k]->dmin);
        runner_do_stars_sort(r, c->progeny[k], flags, cleanup_prog, 0);
Loic Hausammann's avatar
Loic Hausammann committed
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
        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];
For faster browsing, not all history is shown. View entire blame