runner.c 103 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
#include "runner_doiact_stars.h"
113
114
115
116
#undef FUNCTION

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

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

            /* Off we go ! */
            continue;

          } else {

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

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

251
        /* We now have a particle whose smoothing length has converged */
252

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

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

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

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

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

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

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

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

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

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

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
321
  if (timer) TIMER_TOC(timer_dostars_ghost);
322
323
}

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

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

340
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
341

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

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

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

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

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

364
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
365
366
}

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

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

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

  TIMER_TIC;

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

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

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

439
440
441
      /* 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
442

443
      if (part_is_active(p, e)) {
444

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

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

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

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

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

476
  struct engine *e = r->e;
477
  const struct cosmology *cosmo = e->cosmology;
478
  const struct star_formation *starform = e->star_formation;
Folkert Nobels's avatar
Folkert Nobels committed
479
  const struct phys_const *constants = e->physical_constants;
480
481
482
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
483
  const int with_cosmology = (e->policy & engine_policy_cosmology);
484
485
486
  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;
487
488
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
489
490
491
492
493
494
495
496
497
498
499

  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 {
500
501
502
503
504
505
506
507
508
509

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

510
511
        /* Calculate the time step of the current particle for cosmo and no
         * cosmo*/
512
513
514
515
516
517
518
519
520
521
522
523
524
        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);
        }

525
526
527
        if (star_formation_should_convert_to_star(
                e, starform, p, xp, constants, cosmo, hydro_props, us, cooling,
                dt_star, with_cosmology)) {
528
          /* Convert your particle to a star */
529
          struct spart *sp = cell_convert_part_to_spart(e, c, p, xp);
530

531
          /* Copy the properties of the gas particle to the star particle */
532
          star_formation_copy_properties(e, p, xp, sp, starform, constants,
Folkert Nobels's avatar
Folkert Nobels committed
533
                                         cosmo, with_cosmology);
Folkert Nobels's avatar
Folkert Nobels committed
534
        }
535
536
      }
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
537
538
539
540
541
  }

  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
542
543
544
545
546
547
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
548
void runner_do_sort_ascending(struct entry *sort, int N) {
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
595
596
597
598
599
600
601
602

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

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

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
644

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

  struct entry *fingers[8];
660
661
662
  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
663
  float buff[8];
664

665
666
  TIMER_TIC;

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

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

680
681
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
682
  runner_check_sorts_hydro(c, c->hydro.sorted);
683
684

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

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

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

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

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

    /* Loop over the 13 different sort arrays. */
727
    for (int j = 0; j < 13; j++) {
728
729
730
731
732

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

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

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

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

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

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

781
782
783
      } /* Merge. */

      /* Add a sentinel. */
784
785
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
786
787

      /* Mark as sorted. */
788
      atomic_or(&c->hydro.sorted, 1 << j);
789
790
791
792
793
794
795
796

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

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

      /* 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;
        }
811
      }
812
813
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
814
815
    }

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

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

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

850
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
851
  runner_check_sorts_hydro(c, flags);
852
853

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

861
  /* Clear the cell's sort flags. */
862
863
864
  c->hydro.do_sort = 0;
  c->hydro.do_sub_sort = 0;
  c->hydro.requires_sorts = 0;
865

866
867
868
  if (clock) TIMER_TOC(timer_dosort);
}

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