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

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

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

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

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

Pedro Gonnet's avatar
Pedro Gonnet committed
40
/* Local headers. */
41
#include "active.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
42
#include "approx_math.h"
43
#include "atomic.h"
44
#include "cell.h"
45
#include "chemistry.h"
46
#include "const.h"
Stefan Arridge's avatar
Stefan Arridge committed
47
#include "cooling.h"
48
#include "debug.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
49
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
50
#include "engine.h"
51
#include "entropy_floor.h"
52
#include "error.h"
53
54
#include "gravity.h"
#include "hydro.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
55
#include "hydro_properties.h"
56
#include "kick.h"
Loikki's avatar
Loikki committed
57
#include "logger.h"
58
#include "minmax.h"
James Willis's avatar
James Willis committed
59
#include "runner_doiact_vec.h"
60
#include "scheduler.h"
61
#include "sort_part.h"
62
#include "space.h"
63
#include "space_getsid.h"
64
#include "star_formation.h"
65
#include "stars.h"
66
67
#include "task.h"
#include "timers.h"
68
#include "timestep.h"
69
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
70
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
71

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

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

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

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

100
101
102
103
104
105
106
/* Import the limiter loop functions. */
#define FUNCTION limiter
#define FUNCTION_TASK_LOOP TASK_LOOP_LIMITER
#include "runner_doiact.h"
#undef FUNCTION
#undef FUNCTION_TASK_LOOP

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

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

/* Import the stars feedback loop functions. */
#define FUNCTION feedback
Loic Hausammann's avatar
Loic Hausammann committed
119
#include "runner_doiact_stars.h"
120
#undef FUNCTION
121

122
123
int emergency_sort = 0;

124
125
126
127
128
129
130
131
/**
 * @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
132
void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
133

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

  TIMER_TIC;

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

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

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

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

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

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

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

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

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

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

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

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

214
215
            stars_reset_feedback(sp);

216
217
218
219
220
            /* Ok, we are done with this particle */
            continue;
          }

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

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

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
248
            stars_init_spart(sp);
249
250
251
252
253
254
255

            /* Off we go ! */
            continue;

          } else {

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

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

265
        /* We now have a particle whose smoothing length has converged */
266
        stars_reset_feedback(sp);
267

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

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

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

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

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

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

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

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

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

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

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
336
  if (timer) TIMER_TOC(timer_dostars_ghost);
337
338
}

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

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

355
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
356

357
  /* Anything to do here? */
358
  if (!cell_is_active_gravity(c, e)) return;
359

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

366
367
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
368

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

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

379
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
380
381
}

382
383
384
385
386
387
388
389
390
/**
 * @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) {

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

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

  TIMER_TIC;

442
  /* Anything to do here? */
443
  if (!cell_is_active_hydro(c, e)) return;
444

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

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

454
455
456
      /* 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
457

458
      if (part_is_active(p, e)) {
459

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

466
467
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
468
469
470
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

471
472
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
473
          dt_therm = get_timestep(p->time_bin, time_base);
474
        }
475

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

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
486
487
488
/**
 *
 */
489
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
Matthieu Schaller's avatar
Matthieu Schaller committed
490

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

  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++)
514
      if (c->progeny[k] != NULL) runner_do_star_formation(r, c->progeny[k], 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
515
  } else {
516
517
518
519
520
521
522
523

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

524
      /* Only work on active particles */
525
526
      if (part_is_active(p, e)) {

527
528
529
        /* Is this particle star forming? */
        if (star_formation_is_star_forming(p, xp, sf_props, phys_const, cosmo,
                                           hydro_props, us, cooling)) {
530

531
532
533
534
535
536
          /* Time-step size for this particle */
          double dt_star;
          if (with_cosmology) {
            const integertime_t ti_step = get_integer_timestep(p->time_bin);
            const integertime_t ti_begin =
                get_integer_time_begin(ti_current - 1, p->time_bin);
537

538
539
540
541
542
543
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

545
546
547
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
548

549
550
551
552
553
554
555
556
557
558
          /* Are we forming a star particle from this SF rate? */
          if (star_formation_should_convert_to_star(p, xp, sf_props, e,
                                                    dt_star)) {

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

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

560
561
562
563
            message(
                "STAR FORMED!!!! c->ID=%d super->ID=%d c->depth=%d"
                "c->maxdepth=%d",
                c->cellID, c->super->cellID, c->depth, c->maxdepth);
564
565
566
567
568
569
570
571
572
573
574
          }

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

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

        } /* Not Star-forming? */
      }   /* is active? */
    }     /* Loop over particles */
Matthieu Schaller's avatar
Matthieu Schaller committed
575
576
  }

577
  if ((c == c->hydro.super) && (current_stars_count != c->stars.count)) {
578
579
580
581
    message(
        "Emergency sort! c->ID=%d c->depth=%d c->maxdepth=%d c=%p c->super=%p",
        c->cellID, c->depth, c->maxdepth, c, c->hydro.super);

582
583
    cell_clear_stars_sort_flags(c, /*is_super=*/1);
    runner_do_stars_sort(r, c, 0x1FFF, /*cleanup=*/0, /*timer=*/0);
584
  }
585

Matthieu Schaller's avatar
Matthieu Schaller committed
586
587
588
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
589
590
591
592
593
594
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
595
void runner_do_sort_ascending(struct entry *sort, int N) {
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649

  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
650
        }
651
652
653
654
655
656
657
658
659
660
661
662
      } 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
663
    }
664
665
666
  }
}

667
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
668
669
670
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
671
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
672
 */
673
674
675
676
677
678
679
680
681
#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
682
#else
Loic Hausammann's avatar
Loic Hausammann committed
683
684
685
686
#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
687
#endif
Loic Hausammann's avatar
Loic Hausammann committed
688
689
690

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
691

Pedro Gonnet's avatar
Pedro Gonnet committed
692
693
694
695
696
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
697
 * @param flags Cell flag.
698
699
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
700
701
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
702
 */
Loic Hausammann's avatar
Loic Hausammann committed
703
704
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
705
706

  struct entry *fingers[8];
707
708
709
  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
710
  float buff[8];
711

712
713
  TIMER_TIC;

714
  /* We need to do the local sorts plus whatever was requested further up. */
715
  flags |= c->hydro.do_sort;
716
  if (cleanup) {
717
    c->hydro.sorted = 0;
718
  } else {
719
    flags &= ~c->hydro.sorted;
720
  }
721
  if (flags == 0 && !c->hydro.do_sub_sort) return;
722
723

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

727
728
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
729
  runner_check_sorts_hydro(c, c->hydro.sorted);
730
731

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
732
733
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
734
735
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
736
  }
737

738
739
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
740
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
741
#endif
742

743
744
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
745
746
747
    if ((flags & (1 << j)) && c->hydro.sort[j] == NULL) {
      if ((c->hydro.sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                                     (count + 1))) == NULL)
748
749
        error("Failed to allocate sort memory.");
    }
750
751
752
753
754
755
  }

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

    /* Fill in the gaps within the progeny. */
756
    float dx_max_sort = 0.0f;
757
    float dx_max_sort_old = 0.0f;
758
    for (int k = 0; k < 8; k++) {
759
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
760
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
761
        runner_do_hydro_sort(r, c->progeny[k], flags,
762
763
764
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
765
766
767
        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);
768
      }
769
    }
770
771
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
772
773

    /* Loop over the 13 different sort arrays. */
774
    for (int j = 0; j < 13; j++) {
775
776
777
778
779

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

      /* Init the particle index offsets. */
780
      int off[8];
781
782
      off[0] = 0;
      for (int k = 1; k < 8; k++)
783
        if (c->progeny[k - 1] != NULL)
784
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
785
786
787
788
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
789
      int inds[8];
790
      for (int k = 0; k < 8; k++) {
791
        inds[k] = k;
792
793
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
794
795
796
797
798
799
800
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
801
802
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
803
          if (buff[inds[k]] < buff[inds[i]]) {
804
            int temp_i = inds[i];
805
806
807
808
809
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
810
      struct entry *finger = c->hydro.sort[j];
811
      for (int ind = 0; ind < count; ind++) {
812
813
814
815
816
817
818
819
820
821

        /* 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. */
822
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
823
          int temp_i = inds[k - 1];
824
825
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
826
        }
827

828
829
830
      } /* Merge. */

      /* Add a sentinel. */
831
832
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
833
834

      /* Mark as sorted. */
835
      atomic_or(&c->hydro.sorted, 1 << j);
836
837
838
839
840
841
842
843

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

844
    /* Reset the sort distance */
845
    if (c->hydro.sorted == 0) {
846
847
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
848
        error("Have non-NULL xparts in foreign cell");
849
#endif
850
851
852
853
854
855
856
857

      /* 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;
        }
858
      }
859
860
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
861
862
    }

863
    /* Fill the sort array. */
864
    for (int k = 0; k < count; k++) {
865
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
866
      for (int j = 0; j < 13; j++)
867
        if (flags & (1 << j)) {
868
869
870
871
          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];
872
        }
873
    }
874
875

    /* Add the sentinel and sort. */
876
    for (int j = 0; j < 13; j++)
877
      if (flags & (1 << j)) {
878
879
880
881
        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);
882
883
884
      }
  }

885
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
886
  /* Verify the sorting. */
887
  for (int j = 0; j < 13; j++) {
888
    if (!(flags & (1 << j))) continue;
889
    struct entry *finger = c->hydro.sort[j];
890
    for (int k = 1; k < count; k++) {
891
892
893
894
895
      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
896

897
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
898
  runner_check_sorts_hydro(c, flags);
899
900

  /* Make sure the sort flags are consistent (upward). */
Pedro Gonnet's avatar
Pedro Gonnet committed
901
902
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
903
904
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags.");
905
  }
906
#endif
907

908
  /* Clear the cell's sort flags. */
909
910
911
  c->hydro.do_sort = 0;
  c->hydro.do_sub_sort = 0;
  c->hydro.requires_sorts = 0;
912

913
914
915
  if (clock) TIMER_TOC(timer_dosort);
}

Loic Hausammann's avatar
Loic Hausammann committed
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
943
944
945
946
947
/**
 * @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 */
948
  if (flags && !cell_are_spart_drifted(c, r->e)) {
Loic Hausammann's avatar
Loic Hausammann committed
949
    error("Sorting un-drifted cell c->nodeID=%d", c->nodeID);
950
  }
Loic Hausammann's avatar
Loic Hausammann committed
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

#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. */
986
987
988
989
        const int cleanup_prog =
            cleanup && (c->progeny[k]->stars.dx_max_sort_old >
                        space_maxreldx * c->progeny[k]->dmin);
        runner_do_stars_sort(r, c->progeny[k], flags, cleanup_prog, 0);
Loic Hausammann's avatar
Loic Hausammann committed
990
991
992
993
994
995
996
997
998
999
1000
        dx_max_sort = max(dx_max_sort, c->progeny[k]->stars.dx_max_sort);
        dx_max_sort_old =
            max(dx_max_sort_old, c->progeny[k]->stars.dx_max_sort_old);
      }
    }
    c->stars.dx_max_sort = dx_max_sort;
    c->stars.dx_max_sort_old = dx_max_sort_old;

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

For faster browsing, not all history is shown. View entire blame