runner.c 75.5 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"
56
#include "minmax.h"
James Willis's avatar
James Willis committed
57
#include "runner_doiact_vec.h"
58
#include "scheduler.h"
59
#include "sort_part.h"
60
#include "sourceterms.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"
Pedro Gonnet's avatar
Pedro Gonnet committed
67

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

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

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

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

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

99
100
101
/* Import the star loop functions. */
#include "runner_doiact_star.h"

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

  TIMER_TIC;

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

124
    if (count > 0) {
Tom Theuns's avatar
Tom Theuns committed
125

126
127
128
129
130
131
      /* do sourceterms in this cell? */
      const int incell =
          sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
      if (incell == 1) {
        sourceterms_apply(r, sourceterms, c);
      }
Tom Theuns's avatar
Tom Theuns committed
132
133
    }
  }
Tom Theuns's avatar
Tom Theuns committed
134
135
136
137

  if (timer) TIMER_TOC(timer_dosource);
}

138
139
140
141
142
143
144
145
146
147
148
149
150
151
/**
 * @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 ?
 */
void runner_do_star_ghost(struct runner *r, struct cell *c, int timer) {

  struct spart *restrict sparts = c->sparts;
  const struct engine *e = r->e;
  const struct space *s = e->s;
  const struct cosmology *cosmo = e->cosmology;
Loic Hausammann's avatar
Loic Hausammann committed
152
153
154
  const struct stars_props *stars_properties = e->stars_properties;
  const float star_h_max = e->stars_properties->h_max;
  const float eps = e->stars_properties->h_tolerance;
155
  const float star_eta_dim =
Loic Hausammann's avatar
Loic Hausammann committed
156
157
      pow_dimension(e->stars_properties->eta_neighbours);
  const int max_smoothing_iter = e->stars_properties->max_smoothing_iterations;
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
  int redo = 0, count = 0;

  TIMER_TIC;

  /* Anything to do here? */
  if (!cell_is_active_star(c, e)) return;

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_star_ghost(r, c->progeny[k], 0);
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
    if ((sid = (int *)malloc(sizeof(int) * c->scount)) == NULL)
      error("Can't allocate memory for pid.");
    for (int k = 0; k < c->scount; k++)
      if (spart_is_active(&sparts[k], e)) {
        sid[count] = k;
        ++count;
      }

    /* While there are particles that need to be updated... */
    for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
         num_reruns++) {

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

      /* Loop over the remaining active parts in this cell. */
      for (int i = 0; i < count; i++) {

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

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
        if (!spart_is_active(sp, e)) error("Ghost applied to inactive particle");
#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;

        if (sp->wcount == 0.f) { /* No neighbours case */

          /* 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
216
          star_end_density(sp, cosmo);
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253

          /* Compute one step of the Newton-Raphson scheme */
          const float n_sum = sp->wcount * h_old_dim;
          const float n_target = star_eta_dim;
          const float f = n_sum - n_target;
          const float f_prime =
              sp->wcount_dh * h_old_dim +
              hydro_dimension * sp->wcount * h_old_dim_minus_one;

          /* 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 */
          if (sp->h < star_h_max) {

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

            /* Re-initialise everything */
Loic Hausammann's avatar
Loic Hausammann committed
254
            star_init_spart(sp, stars_properties);
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352

            /* Off we go ! */
            continue;

          } else {

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

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
              star_spart_has_no_neighbours(sp, cosmo);
            }
          }
        }

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

        /* As of here, particle force variables will be set. */

        /* Compute variables required for the force loop */
        star_prepare_force(sp, cosmo);

        /* The particle force values are now set.  Do _NOT_
           try to read any particle density variables! */

        /* Prepare the particle for the force loop over neighbours */
        star_reset_acceleration(sp);

      }

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

      /* Re-set the counter for the next loop (potentially). */
      count = redo;
      if (count > 0) {

        /* Climb up the cell hierarchy. */
        for (struct cell *finger = c; finger != NULL; finger = finger->parent) {

          /* Run through this cell's density interactions. */
          for (struct link *l = finger->density; l != NULL; l = l->next) {

#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)
              runner_doself_subset_star_density(r, finger, parts, pid, count);

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

              /* Left or right? */
              if (l->t->ci == finger)
                runner_dopair_subset_star_density(r, finger, parts, pid,
                                                    count, l->t->cj);
              else
                runner_dopair_subset_star_density(r, finger, parts, pid,
                                                    count, l->t->ci);
            }

            /* Otherwise, sub-self interaction? */
            else if (l->t->type == task_type_sub_self)
              runner_dosub_subset_density(r, finger, parts, pid, count, NULL,
                                          -1, 1);

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

              /* Left or right? */
              if (l->t->ci == finger)
                runner_dosub_subset_density(r, finger, parts, pid, count,
                                            l->t->cj, -1, 1);
              else
                runner_dosub_subset_density(r, finger, parts, pid, count,
                                            l->t->ci, -1, 1);
            }
          }
        }
      }
    }

    if (count) {
      error("Smoothing length failed to converge on %i particles.", count);
    }

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

  if (timer) TIMER_TOC(timer_do_star_ghost);
}


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

Matthieu Schaller's avatar
Matthieu Schaller committed
362
363
  struct gpart *restrict gparts = c->gparts;
  const int gcount = c->gcount;
364
365
366
  const struct engine *e = r->e;
  const struct external_potential *potential = e->external_potential;
  const struct phys_const *constants = e->physical_constants;
367
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
368

369
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
370

371
  /* Anything to do here? */
372
  if (!cell_is_active_gravity(c, e)) return;
373

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

380
381
    /* Loop over the gparts in this cell. */
    for (int i = 0; i < gcount; i++) {
382

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

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

393
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
394
395
}

396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
/**
 * @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) {

  struct gpart *restrict gparts = c->gparts;
  const int gcount = c->gcount;
  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
431
/**
432
433
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
434
435
436
437
438
439
440
 *
 * @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) {

441
  const struct engine *e = r->e;
442
443
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
444
445
  const struct cooling_function_data *cooling_func = e->cooling_func;
  const struct phys_const *constants = e->physical_constants;
446
  const struct unit_system *us = e->internal_units;
447
  const double time_base = e->time_base;
448
449
450
451
  const integertime_t ti_current = e->ti_current;
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
  const int count = c->count;
Stefan Arridge's avatar
Stefan Arridge committed
452
453
454

  TIMER_TIC;

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

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

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

467
468
469
      /* 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
470

471
      if (part_is_active(p, e)) {
472

473
474
475
476
477
478
479
480
481
482
        double dt_cool;
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(ti_current + 1, p->time_bin);
          dt_cool =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
        } else {
          dt_cool = get_timestep(p->time_bin, time_base);
        }
483

484
        /* Let's cool ! */
485
        cooling_cool_part(constants, us, cosmo, cooling_func, p, xp, dt_cool);
486
      }
Stefan Arridge's avatar
Stefan Arridge committed
487
488
489
490
491
492
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
493
494
495
496
497
498
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
499
void runner_do_sort_ascending(struct entry *sort, int N) {
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
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

  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
554
        }
555
556
557
558
559
560
561
562
563
564
565
566
      } 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
567
    }
568
569
570
  }
}

Matthieu Schaller's avatar
Matthieu Schaller committed
571
572
573
574
575
576
577
578
/**
 * @brief Recursively checks that the flags are consistent in a cell hierarchy.
 *
 * Debugging function.
 *
 * @param c The #cell to check.
 * @param flags The sorting flags to check.
 */
579
void runner_check_sorts(struct cell *c, int flags) {
Matthieu Schaller's avatar
Matthieu Schaller committed
580
581

#ifdef SWIFT_DEBUG_CHECKS
Pedro Gonnet's avatar
Pedro Gonnet committed
582
  if (flags & ~c->sorted) error("Inconsistent sort flags (downward)!");
583
584
  if (c->split)
    for (int k = 0; k < 8; k++)
585
      if (c->progeny[k] != NULL && c->progeny[k]->count > 0)
586
        runner_check_sorts(c->progeny[k], c->sorted);
Matthieu Schaller's avatar
Matthieu Schaller committed
587
588
589
#else
  error("Calling debugging code without debugging flag activated.");
#endif
590
591
}

Pedro Gonnet's avatar
Pedro Gonnet committed
592
593
594
595
596
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
597
 * @param flags Cell flag.
598
599
 * @param cleanup If true, re-build the sorts for the selected flags instead
 *        of just adding them.
600
601
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
602
 */
603
604
void runner_do_sort(struct runner *r, struct cell *c, int flags, int cleanup,
                    int clock) {
605
606

  struct entry *fingers[8];
607
  const int count = c->count;
608
609
  const struct part *parts = c->parts;
  struct xpart *xparts = c->xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
610
  float buff[8];
611

612
613
  TIMER_TIC;

614
615
  /* We need to do the local sorts plus whatever was requested further up. */
  flags |= c->do_sort;
616
617
618
619
620
  if (cleanup) {
    c->sorted = 0;
  } else {
    flags &= ~c->sorted;
  }
621
  if (flags == 0 && !c->do_sub_sort) return;
622
623

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

627
628
629
630
631
#ifdef SWIFT_DEBUG_CHECKS
  /* Make sure the sort flags are consistent (downward). */
  runner_check_sorts(c, c->sorted);

  /* Make sure the sort flags are consistent (upard). */
Pedro Gonnet's avatar
Pedro Gonnet committed
632
633
634
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
    if (finger->sorted & ~c->sorted) error("Inconsistent sort flags (upward).");
635
  }
636

637
638
  /* Update the sort timer which represents the last time the sorts
     were re-set. */
639
  if (c->sorted == 0) c->ti_sort = r->e->ti_current;
640
#endif
641

642
643
644
645
646
647
648
  /* start by allocating the entry arrays in the requested dimensions. */
  for (int j = 0; j < 13; j++) {
    if ((flags & (1 << j)) && c->sort[j] == NULL) {
      if ((c->sort[j] = (struct entry *)malloc(sizeof(struct entry) *
                                               (count + 1))) == NULL)
        error("Failed to allocate sort memory.");
    }
649
650
651
652
653
654
  }

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

    /* Fill in the gaps within the progeny. */
655
    float dx_max_sort = 0.0f;
656
    float dx_max_sort_old = 0.0f;
657
    for (int k = 0; k < 8; k++) {
658
      if (c->progeny[k] != NULL && c->progeny[k]->count > 0) {
659
660
661
662
663
        /* Only propagate cleanup if the progeny is stale. */
        runner_do_sort(r, c->progeny[k], flags,
                       cleanup && (c->progeny[k]->dx_max_sort >
                                   space_maxreldx * c->progeny[k]->dmin),
                       0);
664
        dx_max_sort = max(dx_max_sort, c->progeny[k]->dx_max_sort);
665
        dx_max_sort_old = max(dx_max_sort_old, c->progeny[k]->dx_max_sort_old);
666
      }
667
    }
668
    c->dx_max_sort = dx_max_sort;
669
    c->dx_max_sort_old = dx_max_sort_old;
670
671

    /* Loop over the 13 different sort arrays. */
672
    for (int j = 0; j < 13; j++) {
673
674
675
676
677

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

      /* Init the particle index offsets. */
678
      int off[8];
679
680
      off[0] = 0;
      for (int k = 1; k < 8; k++)
681
682
683
684
685
686
        if (c->progeny[k - 1] != NULL)
          off[k] = off[k - 1] + c->progeny[k - 1]->count;
        else
          off[k] = off[k - 1];

      /* Init the entries and indices. */
687
      int inds[8];
688
      for (int k = 0; k < 8; k++) {
689
690
        inds[k] = k;
        if (c->progeny[k] != NULL && c->progeny[k]->count > 0) {
691
          fingers[k] = c->progeny[k]->sort[j];
692
693
694
695
696
697
698
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

      /* Sort the buffer. */
699
700
      for (int i = 0; i < 7; i++)
        for (int k = i + 1; k < 8; k++)
701
          if (buff[inds[k]] < buff[inds[i]]) {
702
            int temp_i = inds[i];
703
704
705
706
707
            inds[i] = inds[k];
            inds[k] = temp_i;
          }

      /* For each entry in the new sort list. */
708
      struct entry *finger = c->sort[j];
709
      for (int ind = 0; ind < count; ind++) {
710
711
712
713
714
715
716
717
718
719

        /* 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. */
720
        for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
721
          int temp_i = inds[k - 1];
722
723
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
724
        }
725

726
727
728
      } /* Merge. */

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

      /* Mark as sorted. */
733
      atomic_or(&c->sorted, 1 << j);
734
735
736
737
738
739
740
741

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

742
    /* Reset the sort distance */
743
    if (c->sorted == 0) {
744
745
#ifdef SWIFT_DEBUG_CHECKS
      if (xparts != NULL && c->nodeID != engine_rank)
746
        error("Have non-NULL xparts in foreign cell");
747
#endif
748
749
750
751
752
753
754
755

      /* 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;
        }
756
      }
757
758
      c->dx_max_sort_old = 0.f;
      c->dx_max_sort = 0.f;
759
760
    }

761
    /* Fill the sort array. */
762
    for (int k = 0; k < count; k++) {
763
      const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
764
      for (int j = 0; j < 13; j++)
765
        if (flags & (1 << j)) {
766
767
768
769
          c->sort[j][k].i = k;
          c->sort[j][k].d = px[0] * runner_shift[j][0] +
                            px[1] * runner_shift[j][1] +
                            px[2] * runner_shift[j][2];
770
        }
771
    }
772
773

    /* Add the sentinel and sort. */
774
    for (int j = 0; j < 13; j++)
775
      if (flags & (1 << j)) {
776
777
778
        c->sort[j][count].d = FLT_MAX;
        c->sort[j][count].i = 0;
        runner_do_sort_ascending(c->sort[j], count);
779
        atomic_or(&c->sorted, 1 << j);
780
781
782
      }
  }

783
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
784
  /* Verify the sorting. */
785
  for (int j = 0; j < 13; j++) {
786
    if (!(flags & (1 << j))) continue;
787
    struct entry *finger = c->sort[j];
788
    for (int k = 1; k < count; k++) {
789
790
791
792
793
      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
794

795
796
797
798
  /* Make sure the sort flags are consistent (downward). */
  runner_check_sorts(c, flags);

  /* Make sure the sort flags are consistent (upward). */
Pedro Gonnet's avatar
Pedro Gonnet committed
799
800
801
  for (struct cell *finger = c->parent; finger != NULL;
       finger = finger->parent) {
    if (finger->sorted & ~c->sorted) error("Inconsistent sort flags.");
802
  }
803
#endif
804

805
806
807
808
809
  /* Clear the cell's sort flags. */
  c->do_sort = 0;
  c->do_sub_sort = 0;
  c->requires_sorts = 0;

810
811
812
  if (clock) TIMER_TOC(timer_dosort);
}

813
/**
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
 * @brief Initialize the multipoles before the gravity calculation.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer 1 if the time is to be recorded.
 */
void runner_do_init_grav(struct runner *r, struct cell *c, int timer) {

  const struct engine *e = r->e;

  TIMER_TIC;

#ifdef SWIFT_DEBUG_CHECKS
  if (!(e->policy & engine_policy_self_gravity))
    error("Grav-init task called outside of self-gravity calculation");
#endif

  /* Anything to do here? */
832
  if (!cell_is_active_gravity(c, e)) return;
833
834

  /* Reset the gravity acceleration tensors */
835
  gravity_field_tensors_init(&c->multipole->pot, e->ti_current);
836
837
838
839
840
841
842
843
844
845
846

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) runner_do_init_grav(r, c->progeny[k], 0);
    }
  }

  if (timer) TIMER_TOC(timer_init_grav);
}

847
/**
848
849
850
851
852
 * @brief Intermediate task after the gradient loop that does final operations
 * on the gradient quantities and optionally slope limits the gradients
 *
 * @param r The runner thread.
 * @param c The cell.
853
 * @param timer Are we timing this ?
854
 */
855
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
856

857
#ifdef EXTRA_HYDRO_LOOP
858

859
  struct part *restrict parts = c->parts;
860
  struct xpart *restrict xparts = c->xparts;
861
  const int count = c->count;
862
  const struct engine *e = r->e;
863
  const struct cosmology *cosmo = e->cosmology;
864

865
866
  TIMER_TIC;

867
  /* Anything to do here? */
868
  if (!cell_is_active_hydro(c, e)) return;
869

870
871
872
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
873
      if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k], 0);
874
875
876
877
878
879
880
  } else {

    /* Loop over the parts in this cell. */
    for (int i = 0; i < count; i++) {

      /* Get a direct pointer on the part. */
      struct part *restrict p = &parts[i];
881
      struct xpart *restrict xp = &xparts[i];
882

883
      if (part_is_active(p, e)) {
884

885
        /* Finish the gradient calculation */
886
        hydro_end_gradient(p);
887
888
889
890
891
892
893
894
895
896
897

        /* As of here, particle force variables will be set. */

        /* Compute variables required for the force loop */
        hydro_prepare_force(p, xp, cosmo);

        /* The particle force values are now set.  Do _NOT_
           try to read any particle density variables! */

        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
898
899
900
      }
    }
  }
901

902
903
  if (timer) TIMER_TOC(timer_do_extra_ghost);

904
905
#else
  error("SWIFT was not compiled with the extra hydro loop activated.");
906
#endif
907
}
908

909
/**
910
911
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
912
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
913
 * @param r The runner thread.
914
 * @param c The cell.
915
 * @param timer Are we timing this ?
916
 */
917
void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
918

Matthieu Schaller's avatar
Matthieu Schaller committed
919
920
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
921
  const struct engine *e = r->e;
922
  const struct space *s = e->s;
923
924
  const struct hydro_space *hs = &s->hs;
  const struct cosmology *cosmo = e->cosmology;
lhausamm's avatar
lhausamm committed
925
  const struct chemistry_global_data *chemistry = e->chemistry;
926
  const float hydro_h_max = e->hydro_properties->h_max;
927
928
929
  const float eps = e->hydro_properties->h_tolerance;
  const float hydro_eta_dim =
      pow_dimension(e->hydro_properties->eta_neighbours);
930
  const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
931
  int redo = 0, count = 0;
932

933
934
  TIMER_TIC;

935
  /* Anything to do here? */
936
  if (!cell_is_active_hydro(c, e)) return;
937

938
939
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
940
    for (int k = 0; k < 8; k++)
941
942
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k], 0);
  } else {
943

944
    /* Init the list of active particles that have to be updated. */
945
    int *pid = NULL;
946
    if ((pid = (int *)malloc(sizeof(int) * c->count)) == NULL)
947
      error("Can't allocate memory for pid.");
948
949
950
951
952
    for (int k = 0; k < c->count; k++)
      if (part_is_active(&parts[k], e)) {
        pid[count] = k;
        ++count;
      }
953

954
955
956
    /* While there are particles that need to be updated... */
    for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
         num_reruns++) {
957

958
959
      /* Reset the redo-count. */
      redo = 0;
960

961
      /* Loop over the remaining active parts in this cell. */
962
      for (int i = 0; i < count; i++) {
963

964
        /* Get a direct pointer on the part. */
965
966
        struct part *p = &parts[pid[i]];
        struct xpart *xp = &xparts[pid[i]];
967

968
#ifdef SWIFT_DEBUG_CHECKS
969
        /* Is this part within the timestep? */
970
971
972
        if (!part_is_active(p, e)) error("Ghost applied to inactive particle");
#endif

973
974
975
976
977
        /* Get some useful values */
        const float h_old = p->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;
978
        int has_no_neighbours = 0;
979

980
        if (p->density.wcount == 0.f) { /* No neighbours case */
981

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

985
986
          /* Double h and try again */
          h_new = 2.f * h_old;
987

988
        } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
989

990
          /* Finish the density calculation */
991
          hydro_end_density(p, cosmo);
992
          chemistry_end_density(p, chemistry, cosmo);
993

994
995
996
997
998
999
1000
          /* Compute one step of the Newton-Raphson scheme */
          const float n_sum = p->density.wcount * h_old_dim;
          const float n_target = hydro_eta_dim;
          const float f = n_sum - n_target;
          const float f_prime =
              p->density.wcount_dh * h_old_dim +
              hydro_dimension * p->density.wcount * h_old_dim_minus_one;
For faster browsing, not all history is shown. View entire blame