runner.c 93.2 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 "sourceterms.h"
62
#include "space.h"
63
#include "space_getsid.h"
64
#include "stars.h"
65
66
#include "task.h"
#include "timers.h"
67
#include "timestep.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
68

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

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

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

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

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

100
101
102
103
104
105
106
/* Import the stars density loop functions. */
#define FUNCTION density
#include "runner_doiact_stars.h"
#undef FUNCTION

/* Import the stars feedback loop functions. */
#define FUNCTION feedback
Loic Hausammann's avatar
Loic Hausammann committed
107
#include "runner_doiact_stars.h"
108
109
#undef FUNCTION

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

  TIMER_TIC;

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

132
    if (count > 0) {
Tom Theuns's avatar
Tom Theuns committed
133

134
135
136
137
138
139
      /* do sourceterms in this cell? */
      const int incell =
          sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
      if (incell == 1) {
        sourceterms_apply(r, sourceterms, c);
      }
Tom Theuns's avatar
Tom Theuns committed
140
141
    }
  }
Tom Theuns's avatar
Tom Theuns committed
142
143
144
145

  if (timer) TIMER_TOC(timer_dosource);
}

146
147
148
149
150
151
152
153
/**
 * @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
154
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
155

156
  struct spart *restrict sparts = c->stars.parts;
157
158
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
Loic Hausammann's avatar
Loic Hausammann committed
159
  const struct stars_props *stars_properties = e->stars_properties;
Loic Hausammann's avatar
Loic Hausammann committed
160
  const float stars_h_max = stars_properties->h_max;
161
  const float eps = stars_properties->h_tolerance;
162
  const float stars_eta_dim = pow_dimension(stars_properties->eta_neighbours);
163
164
  const int max_smoothing_iter = stars_properties->max_smoothing_iterations;
  int redo = 0, scount = 0;
165
166
167
168

  TIMER_TIC;

  /* Anything to do here? */
Loic Hausammann's avatar
Loic Hausammann committed
169
  if (!cell_is_active_stars(c, e)) return;
170
171
172
173

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
Loic Hausammann's avatar
Loic Hausammann committed
174
      if (c->progeny[k] != NULL) runner_do_stars_ghost(r, c->progeny[k], 0);
175
176
177
178
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
179
    if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
Loic Hausammann's avatar
Loic Hausammann committed
180
      error("Can't allocate memory for sid.");
181
    for (int k = 0; k < c->stars.count; k++)
182
      if (spart_is_active(&sparts[k], e)) {
183
184
        sid[scount] = k;
        ++scount;
185
186
187
      }

    /* While there are particles that need to be updated... */
188
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
189
190
191
192
193
194
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
195
      for (int i = 0; i < scount; i++) {
196
197
198
199
200
201

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
202
203
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
204
205
206
207
208
209
210
211
212
#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;

213
        if (sp->density.wcount == 0.f) { /* No neighbours case */
214
215
216
217
218
219
220
221
222

          /* 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
223
          stars_end_density(sp, cosmo);
224
225

          /* Compute one step of the Newton-Raphson scheme */
226
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
227
          const float n_target = stars_eta_dim;
228
229
          const float f = n_sum - n_target;
          const float f_prime =
230
231
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252

          /* 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
253
          if (sp->h < stars_h_max) {
254
255
256
257
258
259

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

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
260
            stars_init_spart(sp);
261
262
263
264
265
266
267

            /* Off we go ! */
            continue;

          } else {

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
268
            sp->h = stars_h_max;
269
270
271

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
272
              stars_spart_has_no_neighbours(sp, cosmo);
273
274
275
276
            }
          }
        }

277
        /* We now have a particle whose smoothing length has converged */
278

Loic Hausammann's avatar
Loic Hausammann committed
279
        /* Compute the stellar evolution  */
280
        stars_evolve_spart(sp, stars_properties, cosmo);
281
282
283
284
285
286
      }

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

      /* Re-set the counter for the next loop (potentially). */
287
288
      scount = redo;
      if (scount > 0) {
289
290
291
292
293

        /* 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
294
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
295
296
297
298
299
300
301
302

#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)
303
304
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
305
306
307
308
309
310

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

              /* Left or right? */
              if (l->t->ci == finger)
311
312
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
313
              else
314
315
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
316
317
318
319
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
320
321
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
322
323
324
325
326
327

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

              /* Left or right? */
              if (l->t->ci == finger)
328
329
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
330
              else
331
332
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
333
334
335
336
337
338
            }
          }
        }
      }
    }

339
340
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
341
342
343
344
345
346
    }

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

Matthieu Schaller's avatar
Matthieu Schaller committed
347
  if (timer) TIMER_TOC(timer_dostars_ghost);
348
349
}

Tom Theuns's avatar
Tom Theuns committed
350
351
352
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
353
354
355
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
356
 */
357
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
358

359
360
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
361
362
363
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
364
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
365

366
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
367

368
  /* Anything to do here? */
369
  if (!cell_is_active_gravity(c, e)) return;
370

Tom Theuns's avatar
Tom Theuns committed
371
372
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
373
    for (int k = 0; k < 8; k++)
374
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
375
  } else {
376

377
378
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
379

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

383
      /* Is this part within the time step? */
384
      if (gpart_is_active(gp, e)) {
385
386
        external_gravity_acceleration(time, potential, constants, gp);
      }
387
    }
388
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
389

390
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
391
392
}

393
394
395
396
397
398
399
400
401
/**
 * @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) {

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

438
  const struct engine *e = r->e;
439
440
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
441
442
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
443
  const struct unit_system *us = e->internal_units;
444
  const struct hydro_props *hydro_props = e->hydro_properties;
445
  const double time_base = e->time_base;
446
  const integertime_t ti_current = e->ti_current;
447
448
449
  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
450
451
452

  TIMER_TIC;

453
  /* Anything to do here? */
454
  if (!cell_is_active_hydro(c, e)) return;
455

Stefan Arridge's avatar
Stefan Arridge committed
456
457
458
459
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
460
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
461

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

465
466
467
      /* 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
468

469
      if (part_is_active(p, e)) {
470

471
        double dt_cool, dt_therm;
472
473
474
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
475
476
              get_integer_time_begin(ti_current - 1, p->time_bin);

477
478
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
479
480
481
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

482
483
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
484
          dt_therm = get_timestep(p->time_bin, time_base);
485
        }
486

487
        /* Let's cool ! */
488
489
        cooling_cool_part(constants, us, cosmo, hydro_props, cooling_func, p,
                          xp, dt_cool, dt_therm);
490
      }
Stefan Arridge's avatar
Stefan Arridge committed
491
492
493
494
495
496
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
497
498
499
500
501
502
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

  const struct engine *e = r->e;
503
  const struct cosmology *cosmo = e->cosmology;
504
505
506
  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
507
508
509
510
511
512
513
514
515
516
517

  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 {
518
519
520
521
522
523
524
525
526
527

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

528
529
        const float rho = hydro_get_physical_density(p, cosmo);

530
        // MATTHIEU: Temporary star-formation law
531
532
533
        // Do not use this at home.
        if (rho > 1.5e7 && e->step > 2) {
          message("Removing particle id=%lld rho=%e", p->id, rho);
534
535
536
537
          cell_convert_part_to_gpart(e, c, p, xp);
        }
      }
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
538
539
540
541
542
  }

  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
543
544
545
546
547
548
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
549
void runner_do_sort_ascending(struct entry *sort, int N) {
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603

  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
604
        }
605
606
607
608
609
610
611
612
613
614
615
616
      } 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
617
    }
618
619
620
  }
}

621
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
622
623
624
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
625
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
626
627
628
 *
 * @param c The #cell to check.
 * @param flags The sorting flags to check.
629
 */
630
631
632
633
634
635
636
637
638
#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
639
#else
Loic Hausammann's avatar
Loic Hausammann committed
640
641
642
643
#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
644
#endif
Loic Hausammann's avatar
Loic Hausammann committed
645
646
647

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
648

Pedro Gonnet's avatar
Pedro Gonnet committed
649
650
651
652
653
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
654
 * @param flags Cell flag.
655
656
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
657
658
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
659
 */
Loic Hausammann's avatar
Loic Hausammann committed
660
661
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
662
663

  struct entry *fingers[8];
664
665
666
  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
667
  float buff[8];
668

669
670
  TIMER_TIC;

671
  /* We need to do the local sorts plus whatever was requested further up. */
672
  flags |= c->hydro.do_sort;
673
  if (cleanup) {
674
    c->hydro.sorted = 0;
675
  } else {
676
    flags &= ~c->hydro.sorted;
677
  }
678
  if (flags == 0 && !c->hydro.do_sub_sort) return;
679
680

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

684
685
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
686
  runner_check_sorts_hydro(c, c->hydro.sorted);
687
688

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
689
690
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
691
692
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
693
  }
694

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

700
701
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
702
703
704
    if ((flags & (1 << j)) && c->hydro.sort[j] == NULL) {
      if ((c->hydro.sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                                     (count + 1))) == NULL)
705
706
        error("Failed to allocate sort memory.");
    }
707
708
709
710
711
712
  }

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

    /* Fill in the gaps within the progeny. */
713
    float dx_max_sort = 0.0f;
714
    float dx_max_sort_old = 0.0f;
715
    for (int k = 0; k < 8; k++) {
716
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
717
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
718
        runner_do_hydro_sort(r, c->progeny[k], flags,
719
720
721
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
722
723
724
        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);
725
      }
726
    }
727
728
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
729
730

    /* Loop over the 13 different sort arrays. */
731
    for (int j = 0; j < 13; j++) {
732
733
734
735
736

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

      /* Init the particle index offsets. */
737
      int off[8];
738
739
      off[0] = 0;
      for (int k = 1; k < 8; k++)
740
        if (c->progeny[k - 1] != NULL)
741
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
742
743
744
745
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
746
      int inds[8];
747
      for (int k = 0; k < 8; k++) {
748
        inds[k] = k;
749
750
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
751
752
753
754
755
756
757
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
758
759
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
760
          if (buff[inds[k]] < buff[inds[i]]) {
761
            int temp_i = inds[i];
762
763
764
765
766
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
767
      struct entry *finger = c->hydro.sort[j];
768
      for (int ind = 0; ind < count; ind++) {
769
770
771
772
773
774
775
776
777
778

        /* 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. */
779
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
780
          int temp_i = inds[k - 1];
781
782
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
783
        }
784

785
786
787
      } /* Merge. */

      /* Add a sentinel. */
788
789
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
790
791

      /* Mark as sorted. */
792
      atomic_or(&c->hydro.sorted, 1 << j);
793
794
795
796
797
798
799
800

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

801
    /* Reset the sort distance */
802
    if (c->hydro.sorted == 0) {
803
804
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
805
        error("Have non-NULL xparts in foreign cell");
806
#endif
807
808
809
810
811
812
813
814

      /* 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;
        }
815
      }
816
817
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
818
819
    }

820
    /* Fill the sort array. */
821
    for (int k = 0; k < count; k++) {
822
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
823
      for (int j = 0; j < 13; j++)
824
        if (flags & (1 << j)) {
825
826
827
828
          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];
829
        }
830
    }
831
832

    /* Add the sentinel and sort. */
833
    for (int j = 0; j < 13; j++)
834
      if (flags & (1 << j)) {
835
836
837
838
        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);
839
840
841
      }
  }

842
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
843
  /* Verify the sorting. */
844
  for (int j = 0; j < 13; j++) {
845
    if (!(flags & (1 << j))) continue;
846
    struct entry *finger = c->hydro.sort[j];
847
    for (int k = 1; k < count; k++) {
848
849
850
851
852
      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
853

854
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
855
  runner_check_sorts_hydro(c, flags);
856
857

  /* Make sure the sort flags are consistent (upward). */
Pedro Gonnet's avatar
Pedro Gonnet committed
858
859
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
860
861
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags.");
862
  }
863
#endif
864

865
  /* Clear the cell's sort flags. */
866
867
868
  c->hydro.do_sort = 0;
  c->hydro.do_sub_sort = 0;
  c->hydro.requires_sorts = 0;
869

870
871
872
  if (clock) TIMER_TOC(timer_dosort);
}

Loic Hausammann's avatar
Loic Hausammann committed
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
931
932
933
934
935
936
937
938
939
940
941
942
/**
 * @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
943
                             cleanup && (c->progeny[k]->stars.dx_max_sort_old >
Loic Hausammann's avatar
Loic Hausammann committed
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;
For faster browsing, not all history is shown. View entire blame