runner.c 93 KB
Newer Older
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
5
6
7
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
 *               2016 John A. Regan (john.a.regan@durham.ac.uk)
 *                    Tom Theuns (tom.theuns@durham.ac.uk)
8
 *
9
10
11
12
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
13
 *
14
15
16
17
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
18
 *
19
20
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
21
 *
22
 ******************************************************************************/
Pedro Gonnet's avatar
Pedro Gonnet committed
23

Pedro Gonnet's avatar
Pedro Gonnet committed
24
25
/* Config parameters. */
#include "../config.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
26
27
28
29

/* Some standard headers. */
#include <float.h>
#include <limits.h>
30
#include <stdlib.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
31

32
33
/* MPI headers. */
#ifdef WITH_MPI
34
#include <mpi.h>
35
36
#endif

37
38
39
/* This object's header. */
#include "runner.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
40
/* Local headers. */
41
#include "active.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
42
#include "approx_math.h"
43
#include "atomic.h"
44
#include "cell.h"
45
#include "chemistry.h"
46
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
47
#include "cooling.h"
48
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
49
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
50
#include "engine.h"
51
#include "error.h"
52
53
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
54
#include "hydro_properties.h"
55
#include "kick.h"
Loikki's avatar
Loikki committed
56
#include "logger.h"
57
#include "minmax.h"
James Willis's avatar
James Willis committed
58
#include "runner_doiact_vec.h"
59
#include "scheduler.h"
60
#include "sort_part.h"
61
#include "space.h"
62
#include "space_getsid.h"
63
#include "stars.h"
64
65
#include "task.h"
#include "timers.h"
66
#include "timestep.h"
Folkert Nobels's avatar
Folkert Nobels committed
67
#include "starformation.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
/* Import the stars density loop functions. */
#define FUNCTION density
Loic Hausammann's avatar
Loic Hausammann committed
102
#include "runner_doiact_stars.h"
103
104
105
106
#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
#undef FUNCTION
109
110
111
112
113
114
115
116
117

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

            /* Off we go ! */
            continue;

          } else {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

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

466
  struct engine *e = r->e;
467
  const struct cosmology *cosmo = e->cosmology;
468
  const struct star_formation *starform = e->star_formation;
Folkert Nobels's avatar
Folkert Nobels committed
469
  const struct phys_const *constants = e->physical_constants;
470
471
472
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
473
  const int with_cosmology = (e->policy & engine_policy_cosmology);
Matthieu Schaller's avatar
Matthieu Schaller committed
474
475
476
477
478
479
480
481
482
483
484

  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 {
485
486
487
488
489
490
491
492
493
494

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

495
        //const float rho = hydro_get_physical_density(p, cosmo);
496
        if (star_formation_convert_to_star(starform, p, xp, constants, cosmo) ) {
497
          star_formation_copy_properties(e, c, p, xp, starform, constants, cosmo, with_cosmology);
498
499
500
501
502
503
504
505
        //struct spart *sp =        cell_conert_part_to_spart(c, p, ...);

//
//
//            
//

//        star_formation_copy_properties(p, xp, sp);
Folkert Nobels's avatar
Folkert Nobels committed
506
        }
507
        // MATTHIEU: Temporary star-formation law
508
        // Do not use this at home.
Folkert Nobels's avatar
Folkert Nobels committed
509
510
511
512
        //if (rho > 1.5e7 && e->step > 2) {
        //  message("Removing particle id=%lld rho=%e", p->id, rho);
        //  cell_convert_part_to_gpart(e, c, p, xp);
        //}
513
514
      }
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
515
516
517
518
519
  }

  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
520
521
522
523
524
525
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
526
void runner_do_sort_ascending(struct entry *sort, int N) {
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580

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

598
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
599
600
601
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
602
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
603
 */
604
605
606
607
608
609
610
611
612
#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
613
#else
Loic Hausammann's avatar
Loic Hausammann committed
614
615
616
617
#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
618
#endif
Loic Hausammann's avatar
Loic Hausammann committed
619
620
621

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
622

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

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

643
644
  TIMER_TIC;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

759
760
761
      } /* Merge. */

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

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

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

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

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

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

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

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

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

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

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

844
845
846
  if (clock) TIMER_TOC(timer_dosort);
}

Loic Hausammann's avatar
Loic Hausammann committed
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
/**
 * @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
917
                             cleanup && (c->progeny[k]->stars.dx_max_sort_old >
Loic Hausammann's avatar
Loic Hausammann committed
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
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
        dx_max_sort = max(dx_max_sort, c->progeny[k]->stars.dx_max_sort);
        dx_max_sort_old =
            max(dx_max_sort_old, c->progeny[k]->stars.dx_max_sort_old);
      }
    }
    c->stars.dx_max_sort = dx_max_sort;
    c->stars.dx_max_sort_old = dx_max_sort_old;

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

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

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

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

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

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

        /* Copy the minimum into the new sort array. */
        finger[ind].d = buff[inds[0]];
        finger[ind].i = fingers[inds[0]]->i + off[inds[0]];

        /* Update the buffer. */
        fingers[inds[0]] += 1;
        buff[inds[0]] = fingers[inds[0]]->d;

        /* Find the smallest entry. */
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
          int temp_i = inds[k - 1];
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
        }

      } /* Merge. */

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

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

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

    /* Reset the sort distance */
    if (c->stars.sorted == 0) {
For faster browsing, not all history is shown. View entire blame