runner.c 109 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 "memuse.h"
59
#include "minmax.h"
James Willis's avatar
James Willis committed
60
#include "runner_doiact_vec.h"
61
#include "scheduler.h"
62
#include "sort_part.h"
63
#include "space.h"
64
#include "space_getsid.h"
65
#include "star_formation.h"
66
#include "star_formation_iact.h"
67
#include "stars.h"
68
69
#include "task.h"
#include "timers.h"
70
#include "timestep.h"
71
#include "timestep_limiter.h"
Folkert Nobels's avatar
Folkert Nobels committed
72
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
73

74
75
76
77
#define TASK_LOOP_DENSITY 0
#define TASK_LOOP_GRADIENT 1
#define TASK_LOOP_FORCE 2
#define TASK_LOOP_LIMITER 3
Loic Hausammann's avatar
Loic Hausammann committed
78
#define TASK_LOOP_FEEDBACK 4
79

80
/* Import the density loop functions. */
81
#define FUNCTION density
82
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
83
#include "runner_doiact.h"
84
85
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
86

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

96
/* Import the force loop functions. */
97
#define FUNCTION force
98
#define FUNCTION_TASK_LOOP TASK_LOOP_FORCE
99
#include "runner_doiact.h"
100
101
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
102

103
104
105
106
107
108
109
/* 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

110
/* Import the gravity loop functions. */
111
#include "runner_doiact_grav.h"
112

113
114
/* Import the stars density loop functions. */
#define FUNCTION density
Loic Hausammann's avatar
Loic Hausammann committed
115
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
Loic Hausammann's avatar
Loic Hausammann committed
116
#include "runner_doiact_stars.h"
Loic Hausammann's avatar
Loic Hausammann committed
117
#undef FUNCTION_TASK_LOOP
118
119
120
121
#undef FUNCTION

/* Import the stars feedback loop functions. */
#define FUNCTION feedback
Loic Hausammann's avatar
Loic Hausammann committed
122
#define FUNCTION_TASK_LOOP TASK_LOOP_FEEDBACK
Loic Hausammann's avatar
Loic Hausammann committed
123
#include "runner_doiact_stars.h"
Loic Hausammann's avatar
Loic Hausammann committed
124
#undef FUNCTION_TASK_LOOP
125
#undef FUNCTION
126
127
128
129
130
131
132
133
134

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

137
  struct spart *restrict sparts = c->stars.parts;
138
139
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
140
141
142
143
144
145
  const float stars_h_max = e->hydro_properties->h_max;
  const float stars_h_min = e->hydro_properties->h_min;
  const float eps = e->stars_properties->h_tolerance;
  const float stars_eta_dim =
      pow_dimension(e->stars_properties->eta_neighbours);
  const int max_smoothing_iter = e->stars_properties->max_smoothing_iterations;
146
  int redo = 0, scount = 0;
147
148
149
150

  TIMER_TIC;

  /* Anything to do here? */
151
  if (c->stars.count == 0) return;
Loic Hausammann's avatar
Loic Hausammann committed
152
  if (!cell_is_active_stars(c, e)) return;
153
154
155
156

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
Loic Hausammann's avatar
Loic Hausammann committed
157
      if (c->progeny[k] != NULL) runner_do_stars_ghost(r, c->progeny[k], 0);
158
159
160
161
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
162
163
164
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
165
    if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
Loic Hausammann's avatar
Loic Hausammann committed
166
      error("Can't allocate memory for sid.");
167
    if ((h_0 = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
168
      error("Can't allocate memory for h_0.");
169
    if ((left = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
170
      error("Can't allocate memory for left.");
171
    if ((right = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
172
      error("Can't allocate memory for right.");
173
    for (int k = 0; k < c->stars.count; k++)
174
      if (spart_is_active(&sparts[k], e)) {
175
        sid[scount] = k;
176
177
178
        h_0[scount] = sparts[k].h;
        left[scount] = 0.f;
        right[scount] = stars_h_max;
179
        ++scount;
180
181
182
      }

    /* While there are particles that need to be updated... */
183
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
184
185
186
187
188
189
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
190
      for (int i = 0; i < scount; i++) {
191
192
193
194
195
196

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
197
198
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
199
200
201
#endif

        /* Get some useful values */
202
        const float h_init = h_0[i];
203
204
205
        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);
206

207
208
209
        float h_new;
        int has_no_neighbours = 0;

210
        if (sp->density.wcount == 0.f) { /* No neighbours case */
211
212
213
214
215
216

          /* Flag that there were no neighbours */
          has_no_neighbours = 1;

          /* Double h and try again */
          h_new = 2.f * h_old;
217

218
219
220
        } else {

          /* Finish the density calculation */
Loic Hausammann's avatar
Loic Hausammann committed
221
          stars_end_density(sp, cosmo);
222
223

          /* Compute one step of the Newton-Raphson scheme */
224
          const float n_sum = sp->density.wcount * h_old_dim;
Loic Hausammann's avatar
Loic Hausammann committed
225
          const float n_target = stars_eta_dim;
226
227
          const float f = n_sum - n_target;
          const float f_prime =
228
229
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
230

231
          /* Improve the bisection bounds */
232
233
234
235
          if (n_sum < n_target)
            left[i] = max(left[i], h_old);
          else if (n_sum > n_target)
            right[i] = min(right[i], h_old);
236
237
238
239
240
241
242

#ifdef SWIFT_DEBUG_CHECKS
          /* Check the validity of the left and right bounds */
          if (left[i] > right[i])
            error("Invalid left (%e) and right (%e)", left[i], right[i]);
#endif

243
          /* Skip if h is already h_max and we don't have enough neighbours */
244
245
246
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {
247

248
249
            stars_reset_feedback(sp);

250
251
252
253
254
            /* Ok, we are done with this particle */
            continue;
          }

          /* Normal case: Use Newton-Raphson to get a better value of h */
255

256
257
          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);
258
259
260
261
262
263
264
265
266
267
268
269

          /* Be verbose about the particles that struggle to converge */
          if (num_reruns > max_smoothing_iter - 10) {

            message(
                "Smoothing length convergence problem: iter=%d p->id=%lld "
                "h_init=%12.8e h_old=%12.8e h_new=%12.8e f=%f f_prime=%f "
                "n_sum=%12.8e n_target=%12.8e left=%12.8e right=%12.8e",
                num_reruns, sp->id, h_init, h_old, h_new, f, f_prime, n_sum,
                n_target, left[i], right[i]);
          }

270
271
272
273
274
275
276
277
278
#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);
279
280
281
282

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
283
284
285
286
287
288
        }

        /* Check whether the particle has an inappropriate smoothing length */
        if (fabsf(h_new - h_old) > eps * h_old) {

          /* Ok, correct then */
289
290
291
292
293
294
295
296
297
298
299
300
301
302

          /* Case where we have been oscillating around the solution */
          if ((h_new == left[i] && h_old == right[i]) ||
              (h_old == left[i] && h_new == right[i])) {

            /* Bissect the remaining interval */
            sp->h = pow_inv_dimension(
                0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));

          } else {

            /* Normal case */
            sp->h = h_new;
          }
303
304

          /* If below the absolute maximum, try again */
305
          if (sp->h < stars_h_max && sp->h > stars_h_min) {
306
307
308

            /* Flag for another round of fun */
            sid[redo] = sid[i];
309
310
311
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
312
313
314
            redo += 1;

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
315
            stars_init_spart(sp);
316
317
318
319

            /* Off we go ! */
            continue;

320
321
322
323
324
325
          } else if (sp->h <= stars_h_min) {

            /* Ok, this particle is a lost cause... */
            sp->h = stars_h_min;

          } else if (sp->h >= stars_h_max) {
326
327

            /* Ok, this particle is a lost cause... */
Loic Hausammann's avatar
Loic Hausammann committed
328
            sp->h = stars_h_max;
329
330
331

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
Loic Hausammann's avatar
Loic Hausammann committed
332
              stars_spart_has_no_neighbours(sp, cosmo);
333
            }
334
335
336
337
338

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
339
340
341
          }
        }

342
        /* We now have a particle whose smoothing length has converged */
343
        stars_reset_feedback(sp);
344

Loic Hausammann's avatar
Loic Hausammann committed
345
        /* Compute the stellar evolution  */
346
        stars_evolve_spart(sp, e->stars_properties, cosmo);
347
348
349
350
351
352
      }

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

      /* Re-set the counter for the next loop (potentially). */
353
354
      scount = redo;
      if (scount > 0) {
355
356
357
358
359

        /* 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
360
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
361
362
363
364
365
366
367
368

#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)
369
370
              runner_doself_subset_branch_stars_density(r, finger, sparts, sid,
                                                        scount);
371
372
373
374
375
376

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

              /* Left or right? */
              if (l->t->ci == finger)
377
378
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->cj);
379
              else
380
381
                runner_dopair_subset_branch_stars_density(
                    r, finger, sparts, sid, scount, l->t->ci);
382
383
384
385
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
386
387
              runner_dosub_subset_stars_density(r, finger, sparts, sid, scount,
                                                NULL, -1, 1);
388
389
390
391
392
393

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

              /* Left or right? */
              if (l->t->ci == finger)
394
395
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->cj, -1, 1);
396
              else
397
398
                runner_dosub_subset_stars_density(r, finger, sparts, sid,
                                                  scount, l->t->ci, -1, 1);
399
400
401
402
403
404
            }
          }
        }
      }
    }

405
406
    if (scount) {
      error("Smoothing length failed to converge on %i particles.", scount);
407
408
409
    }

    /* Be clean */
410
411
    free(left);
    free(right);
412
    free(sid);
413
    free(h_0);
414
415
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
416
  if (timer) TIMER_TOC(timer_dostars_ghost);
417
418
}

Tom Theuns's avatar
Tom Theuns committed
419
420
421
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
422
423
424
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
425
 */
426
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
427

428
429
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
430
431
432
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
433
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
434

435
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
436

437
  /* Anything to do here? */
438
  if (!cell_is_active_gravity(c, e)) return;
439

Tom Theuns's avatar
Tom Theuns committed
440
441
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
442
    for (int k = 0; k < 8; k++)
443
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
444
  } else {
445

446
447
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
448

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

452
      /* Is this part within the time step? */
453
      if (gpart_is_active(gp, e)) {
454
455
        external_gravity_acceleration(time, potential, constants, gp);
      }
456
    }
457
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
458

459
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
460
461
}

462
463
464
465
466
467
468
469
470
/**
 * @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) {

471
472
  struct gpart *restrict gparts = c->grav.parts;
  const int gcount = c->grav.count;
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
  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
497
/**
498
499
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
500
501
502
503
504
505
506
 *
 * @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) {

507
  const struct engine *e = r->e;
508
509
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
510
511
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
512
  const struct unit_system *us = e->internal_units;
513
  const struct hydro_props *hydro_props = e->hydro_properties;
514
  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
515
  const double time_base = e->time_base;
516
  const integertime_t ti_current = e->ti_current;
517
518
519
  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
520
521
522

  TIMER_TIC;

523
  /* Anything to do here? */
524
  if (!cell_is_active_hydro(c, e)) return;
525

Stefan Arridge's avatar
Stefan Arridge committed
526
527
528
529
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
530
  } else {
Stefan Arridge's avatar
Stefan Arridge committed
531

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

535
536
537
      /* 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
538

539
      if (part_is_active(p, e)) {
540

541
        double dt_cool, dt_therm;
542
543
544
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
545
546
              get_integer_time_begin(ti_current - 1, p->time_bin);

547
548
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
549
550
551
          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
                                                     ti_begin + ti_step);

552
553
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
554
          dt_therm = get_timestep(p->time_bin, time_base);
555
        }
556

557
        /* Let's cool ! */
558
559
560
        cooling_cool_part(constants, us, cosmo, hydro_props,
                          entropy_floor_props, cooling_func, p, xp, dt_cool,
                          dt_therm);
561
      }
Stefan Arridge's avatar
Stefan Arridge committed
562
563
564
565
566
567
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Matthieu Schaller's avatar
Matthieu Schaller committed
568
569
570
571
572
/**
 *
 */
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {

573
  struct engine *e = r->e;
574
  const struct cosmology *cosmo = e->cosmology;
575
576
  const struct star_formation *sf_props = e->star_formation;
  const struct phys_const *phys_const = e->physical_constants;
577
578
579
  const int count = c->hydro.count;
  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
580
  const int with_cosmology = (e->policy & engine_policy_cosmology);
581
  const int with_feedback = (e->policy & engine_policy_feedback);
582
583
584
  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;
585
  const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
586
587
  const double time_base = e->time_base;
  const integertime_t ti_current = e->ti_current;
588
  const int current_stars_count = c->stars.count;
Matthieu Schaller's avatar
Matthieu Schaller committed
589
590
591
592
593
594
595
596
597
598
599

  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 {
600
601
602
603
604
605
606
607

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

608
      /* Only work on active particles */
609
610
      if (part_is_active(p, e)) {

611
612
        /* Is this particle star forming? */
        if (star_formation_is_star_forming(p, xp, sf_props, phys_const, cosmo,
Folkert Nobels's avatar
Folkert Nobels committed
613
614
                                           hydro_props, us, cooling,
                                           entropy_floor)) {
615

616
617
618
619
620
621
          /* 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);
622

623
624
625
626
627
628
            dt_star =
                cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);

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

630
631
632
          /* Compute the SF rate of the particle */
          star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
                                     dt_star);
633

634
635
636
637
638
639
640
          /* 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);

641
642
643
644
645
646
647
            /* Did we get a star? (Or did we run out of spare ones?) */
            if (sp != NULL) {

              /* Copy the properties of the gas particle to the star particle */
              star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
                                             with_cosmology);
            }
648
649
650
651
652
653
654
655
656
657
658
          }

        } 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
659
660
  }

661
662
  /* If we formed any stars, the star sorts are now invalid. We need to
   * re-compute them. */
663
664
  if (with_feedback && (c == c->hydro.super) &&
      (current_stars_count != c->stars.count)) {
665
666

    cell_clear_stars_sort_flags(c);
667
    runner_do_stars_sort(r, c, 0x1FFF, /*cleanup=*/0, /*timer=*/0);
668
  }
669

Matthieu Schaller's avatar
Matthieu Schaller committed
670
671
672
  if (timer) TIMER_TOC(timer_do_star_formation);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
673
674
675
676
677
678
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
679
void runner_do_sort_ascending(struct entry *sort, int N) {
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733

  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
734
        }
735
736
737
738
739
740
741
742
743
744
745
746
      } 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
747
    }
748
749
750
  }
}

751
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
752
753
754
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
755
 * Debugging function. Exists in two flavours: hydro & stars.
Matthieu Schaller's avatar
Matthieu Schaller committed
756
 */
757
758
759
760
761
762
763
764
765
#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
766
#else
Loic Hausammann's avatar
Loic Hausammann committed
767
768
769
770
#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
771
#endif
Loic Hausammann's avatar
Loic Hausammann committed
772
773
774

RUNNER_CHECK_SORTS(hydro)
RUNNER_CHECK_SORTS(stars)
775

Pedro Gonnet's avatar
Pedro Gonnet committed
776
777
778
779
780
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
781
 * @param flags Cell flag.
782
783
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
784
785
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
786
 */
Loic Hausammann's avatar
Loic Hausammann committed
787
788
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
                          int cleanup, int clock) {
789
790

  struct entry *fingers[8];
791
792
793
  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
794
  float buff[8];
795

796
797
  TIMER_TIC;

798
  /* We need to do the local sorts plus whatever was requested further up. */
799
  flags |= c->hydro.do_sort;
800
  if (cleanup) {
801
    c->hydro.sorted = 0;
802
  } else {
803
    flags &= ~c->hydro.sorted;
804
  }
805
  if (flags == 0 && !c->hydro.do_sub_sort) return;
806
807

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

811
812
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
813
  runner_check_sorts_hydro(c, c->hydro.sorted);
814
815

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
816
817
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
818
819
    if (finger->hydro.sorted & ~c->hydro.sorted)
      error("Inconsistent sort flags (upward).");
820
  }
821

822
823
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
824
  if (c->hydro.sorted == 0) c->hydro.ti_sort = r->e->ti_current;
825
#endif
826

827
828
  /* Allocate memory for sorting. */
  cell_malloc_hydro_sorts(c, flags);
829
830
831
832
833

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

    /* Fill in the gaps within the progeny. */
834
    float dx_max_sort = 0.0f;
835
    float dx_max_sort_old = 0.0f;
836
    for (int k = 0; k < 8; k++) {
837
      if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
838
        /* Only propagate cleanup if the progeny is stale. */
Loic Hausammann's avatar
Loic Hausammann committed
839
        runner_do_hydro_sort(r, c->progeny[k], flags,
840
841
842
                             cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
                                         space_maxreldx * c->progeny[k]->dmin),
                             0);
843
844
845
        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);
846
      }
847
    }
848
849
    c->hydro.dx_max_sort = dx_max_sort;
    c->hydro.dx_max_sort_old = dx_max_sort_old;
850
851

    /* Loop over the 13 different sort arrays. */
852
    for (int j = 0; j < 13; j++) {
853
854
855
856
857

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

      /* Init the particle index offsets. */
858
      int off[8];
859
860
      off[0] = 0;
      for (int k = 1; k < 8; k++)
861
        if (c->progeny[k - 1] != NULL)
862
          off[k] = off[k - 1] + c->progeny[k - 1]->hydro.count;
863
864
865
866
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
867
      int inds[8];
868
      for (int k = 0; k < 8; k++) {
869
        inds[k] = k;
870
871
        if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
          fingers[k] = c->progeny[k]->hydro.sort[j];
872
873
874
875
876
877
878
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
879
880
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
881
          if (buff[inds[k]] < buff[inds[i]]) {
882
            int temp_i = inds[i];
883
884
885
886
887
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
888
      struct entry *finger = c->hydro.sort[j];
889
      for (int ind = 0; ind < count; ind++) {
890
891
892
893
894
895
896
897
898
899

        /* 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. */
900
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
901
          int temp_i = inds[k - 1];
902
903
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
904
        }
905

906
907
908
      } /* Merge. */

      /* Add a sentinel. */
909
910
      c->hydro.sort[j][count].d = FLT_MAX;
      c->hydro.sort[j][count].i = 0;
911
912

      /* Mark as sorted. */
913
      atomic_or(&c->hydro.sorted, 1 << j);
914
915
916
917
918
919
920
921

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

922
    /* Reset the sort distance */
923
    if (c->hydro.sorted == 0) {
924
925
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
926
        error("Have non-NULL xparts in foreign cell");
927
#endif
928
929
930
931
932
933
934
935

      /* 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;
        }
936
      }
937
938
      c->hydro.dx_max_sort_old = 0.f;
      c->hydro.dx_max_sort = 0.f;
939
940
    }

941
    /* Fill the sort array. */
942
    for (int k = 0; k < count; k++) {
943
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
944
      for (int j = 0; j < 13; j++)
945
        if (flags & (1 << j)) {
946
947
948
949
          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];
950
        }
951
    }
952
953

    /* Add the sentinel and sort. */
954
    for (int j = 0; j < 13; j++)
955
      if (flags & (1 << j)) {
956
957
958
959
        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);
960
961
962
      }
  }

963
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
964
  /* Verify the sorting. */
965
  for (int j = 0; j < 13; j++) {
966
    if (!(flags & (1 << j))) continue;
967
    struct entry *finger = c->hydro.sort[j];
968
    for (int k = 1; k < count; k++) {
969
970
971
972
973
      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
974

975
  /* Make sure the sort flags are consistent (downward). */
Loic Hausammann's avatar
Loic Hausammann committed
976
  runner_check_sorts_hydro(c, flags);