runner.c 31.8 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. */
Matthieu Schaller's avatar
Matthieu Schaller committed
41
#include "approx_math.h"
42
#include "atomic.h"
43
#include "cell.h"
44
#include "const.h"
45
#include "debug.h"
46
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
47
#include "engine.h"
48
#include "error.h"
49
50
#include "gravity.h"
#include "hydro.h"
51
#include "hydro_properties.h"
52
#include "kick.h"
53
#include "minmax.h"
54
55
56
57
#include "scheduler.h"
#include "space.h"
#include "task.h"
#include "timers.h"
58
#include "timestep.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
59

60
/* Orientation of the cell pairs */
61
62
63
64
65
66
67
68
69
70
71
72
73
const double runner_shift[13][3] = {
    {5.773502691896258e-01, 5.773502691896258e-01, 5.773502691896258e-01},
    {7.071067811865475e-01, 7.071067811865475e-01, 0.0},
    {5.773502691896258e-01, 5.773502691896258e-01, -5.773502691896258e-01},
    {7.071067811865475e-01, 0.0, 7.071067811865475e-01},
    {1.0, 0.0, 0.0},
    {7.071067811865475e-01, 0.0, -7.071067811865475e-01},
    {5.773502691896258e-01, -5.773502691896258e-01, 5.773502691896258e-01},
    {7.071067811865475e-01, -7.071067811865475e-01, 0.0},
    {5.773502691896258e-01, -5.773502691896258e-01, -5.773502691896258e-01},
    {0.0, 7.071067811865475e-01, 7.071067811865475e-01},
    {0.0, 1.0, 0.0},
    {0.0, 7.071067811865475e-01, -7.071067811865475e-01},
Pedro Gonnet's avatar
Pedro Gonnet committed
74
75
    {0.0, 0.0, 1.0},
};
76
77

/* Does the axis need flipping ? */
78
79
const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
                              0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
Pedro Gonnet's avatar
Pedro Gonnet committed
80

81
/* Import the density loop functions. */
82
83
84
#define FUNCTION density
#include "runner_doiact.h"

85
/* Import the force loop functions. */
86
87
88
89
#undef FUNCTION
#define FUNCTION force
#include "runner_doiact.h"

Tom Theuns's avatar
Tom Theuns committed
90
91
92
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
93
94
95
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
96
 */
97
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
98

99
100
  struct gpart *restrict gparts = c->gparts;
  const int gcount = c->gcount;
101
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
102
  const struct external_potential *potential = r->e->external_potential;
103
  const struct phys_const *constants = r->e->physical_constants;
Matthieu Schaller's avatar
Matthieu Schaller committed
104

105
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
106
107
108

  /* Recurse? */
  if (c->split) {
109
    for (int k = 0; k < 8; k++)
110
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
Tom Theuns's avatar
Tom Theuns committed
111
112
113
114
    return;
  }

#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
115
  OUT;
Tom Theuns's avatar
Tom Theuns committed
116
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
117

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

121
    /* Get a direct pointer on the part. */
122
    struct gpart *const g = &gparts[i];
123
124
125

    /* Is this part within the time step? */
    if (g->ti_end <= ti_current) {
Matthieu Schaller's avatar
Matthieu Schaller committed
126

127
      external_gravity(potential, constants, g);
128
    }
129
  }
130
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
131
132
}

Pedro Gonnet's avatar
Pedro Gonnet committed
133
134
135
136
137
138
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
139
void runner_do_sort_ascending(struct entry *sort, int N) {
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
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

  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
194
        }
195
196
197
198
199
200
201
202
203
204
205
206
      } 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
207
    }
208
209
210
  }
}

Pedro Gonnet's avatar
Pedro Gonnet committed
211
212
213
214
215
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
216
 * @param flags Cell flag.
217
218
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
219
 */
220
void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
221
222
223
224
225
226
227

  struct entry *finger;
  struct entry *fingers[8];
  struct part *parts = c->parts;
  struct entry *sort;
  int j, k, count = c->count;
  int i, ind, off[8], inds[8], temp_i, missing;
228
229
  float buff[8];
  double px[3];
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253

  TIMER_TIC

  /* Clean-up the flags, i.e. filter out what's already been sorted. */
  flags &= ~c->sorted;
  if (flags == 0) return;

  /* start by allocating the entry arrays. */
  if (c->sort == NULL || c->sortsize < count) {
    if (c->sort != NULL) free(c->sort);
    c->sortsize = count * 1.1;
    if ((c->sort = (struct entry *)malloc(sizeof(struct entry) *
                                          (c->sortsize + 1) * 13)) == NULL)
      error("Failed to allocate sort memory.");
  }
  sort = c->sort;

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

    /* Fill in the gaps within the progeny. */
    for (k = 0; k < 8; k++) {
      if (c->progeny[k] == NULL) continue;
      missing = flags & ~c->progeny[k]->sorted;
254
      if (missing) runner_do_sort(r, c->progeny[k], missing, 0);
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
    }

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

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

      /* Init the particle index offsets. */
      for (off[0] = 0, k = 1; k < 8; k++)
        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. */
      for (k = 0; k < 8; k++) {
        inds[k] = k;
        if (c->progeny[k] != NULL && c->progeny[k]->count > 0) {
          fingers[k] = &c->progeny[k]->sort[j * (c->progeny[k]->count + 1)];
          buff[k] = fingers[k]->d;
          off[k] = off[k];
        } else
          buff[k] = FLT_MAX;
      }

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

      /* For each entry in the new sort list. */
      finger = &sort[j * (count + 1)];
      for (ind = 0; ind < count; ind++) {

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

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

        /* Find the smallest entry. */
        for (k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
          temp_i = inds[k - 1];
          inds[k - 1] = inds[k];
          inds[k] = temp_i;
Pedro Gonnet's avatar
Pedro Gonnet committed
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
      } /* Merge. */

      /* Add a sentinel. */
      sort[j * (count + 1) + count].d = FLT_MAX;
      sort[j * (count + 1) + count].i = 0;

      /* Mark as sorted. */
      c->sorted |= (1 << j);

    } /* loop over sort arrays. */

  } /* progeny? */

  /* Otherwise, just sort. */
  else {

    /* Fill the sort array. */
    for (k = 0; k < count; k++) {
      px[0] = parts[k].x[0];
      px[1] = parts[k].x[1];
      px[2] = parts[k].x[2];
      for (j = 0; j < 13; j++)
        if (flags & (1 << j)) {
          sort[j * (count + 1) + k].i = k;
333
334
335
          sort[j * (count + 1) + k].d = px[0] * runner_shift[j][0] +
                                        px[1] * runner_shift[j][1] +
                                        px[2] * runner_shift[j][2];
336
        }
337
    }
338
339
340
341
342
343

    /* Add the sentinel and sort. */
    for (j = 0; j < 13; j++)
      if (flags & (1 << j)) {
        sort[j * (count + 1) + count].d = FLT_MAX;
        sort[j * (count + 1) + count].i = 0;
344
        runner_do_sort_ascending(&sort[j * (count + 1)], count);
345
346
347
348
        c->sorted |= (1 << j);
      }
  }

349
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
350
  /* Verify the sorting. */
351
352
353
354
355
356
357
  for (j = 0; j < 13; j++) {
    if (!(flags & (1 << j))) continue;
    finger = &sort[j * (count + 1)];
    for (k = 1; k < count; k++) {
      if (finger[k].d < finger[k - 1].d)
        error("Sorting failed, ascending array.");
      if (finger[k].i >= count) error("Sorting failed, indices borked.");
358
359
    }
  }
360
#endif
361
362
363
364

  if (clock) TIMER_TOC(timer_dosort);
}

365
366
367
368
369
/**
 * @brief Initialize the particles before the density calculation
 *
 * @param r The runner thread.
 * @param c The cell.
Matthieu Schaller's avatar
Matthieu Schaller committed
370
 * @param timer 1 if the time is to be recorded.
371
 */
372
void runner_do_init(struct runner *r, struct cell *c, int timer) {
373

374
375
  struct part *const parts = c->parts;
  struct gpart *const gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
376
  const int count = c->count;
377
  const int gcount = c->gcount;
378
  const int ti_current = r->e->ti_current;
379
380

  TIMER_TIC;
381

382
383
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
384
    for (int k = 0; k < 8; k++)
385
      if (c->progeny[k] != NULL) runner_do_init(r, c->progeny[k], 0);
386
    return;
387
388
  } else {

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

392
      /* Get a direct pointer on the part. */
393
      struct part *const p = &parts[i];
394

395
      if (p->ti_end <= ti_current) {
Matthieu Schaller's avatar
Matthieu Schaller committed
396

397
398
        /* Get ready for a density calculation */
        hydro_init_part(p);
399
      }
400
    }
401
402
403
404
405
406
407
408
409
410
411
412
413

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

      /* Get a direct pointer on the part. */
      struct gpart *const gp = &gparts[i];

      if (gp->ti_end <= ti_current) {

        /* Get ready for a density calculation */
        gravity_init_part(gp);
      }
    }
414
  }
415

Peter W. Draper's avatar
Peter W. Draper committed
416
  if (timer) TIMER_TOC(timer_init);
417
418
}

419
420
421
/**
 * @brief Intermediate task between density and force
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
422
 * @param r The runner thread.
423
 * @param c The cell.
424
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
425
void runner_do_ghost(struct runner *r, struct cell *c) {
426
427

  struct part *p, *parts = c->parts;
428
  struct xpart *xp, *xparts = c->xparts;
429
  struct cell *finger;
Matthieu Schaller's avatar
Matthieu Schaller committed
430
  int redo, count = c->count;
431
  int *pid;
432
  float h_corr;
433
434
  const int ti_current = r->e->ti_current;
  const double timeBase = r->e->timeBase;
435
436
437
438
439
440
441
  const float target_wcount = r->e->hydro_properties->target_neighbours;
  const float max_wcount =
      target_wcount + r->e->hydro_properties->delta_neighbours;
  const float min_wcount =
      target_wcount - r->e->hydro_properties->delta_neighbours;
  const int max_smoothing_iter =
      r->e->hydro_properties->max_smoothing_iterations;
442

443
444
  TIMER_TIC;

445
446
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
447
    for (int k = 0; k < 8; k++)
Matthieu Schaller's avatar
Matthieu Schaller committed
448
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k]);
449
450
451
452
453
454
    return;
  }

  /* Init the IDs that have to be updated. */
  if ((pid = (int *)alloca(sizeof(int) * count)) == NULL)
    error("Call to alloca failed.");
Matthieu Schaller's avatar
Matthieu Schaller committed
455
  for (int k = 0; k < count; k++) pid[k] = k;
456
457

  /* While there are particles that need to be updated... */
458
  for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
Matthieu Schaller's avatar
Matthieu Schaller committed
459
       num_reruns++) {
460
461
462

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

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

      /* Get a direct pointer on the part. */
      p = &parts[pid[i]];
469
      xp = &xparts[pid[i]];
470
471

      /* Is this part within the timestep? */
472
      if (p->ti_end <= ti_current) {
473

474
        /* Finish the density calculation */
475
        hydro_end_density(p, ti_current);
476
477

        /* If no derivative, double the smoothing length. */
478
        if (p->density.wcount_dh == 0.0f) h_corr = p->h;
479
480
481

        /* Otherwise, compute the smoothing length update (Newton step). */
        else {
482
          h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
483
484

          /* Truncate to the range [ -p->h/2 , p->h ]. */
485
486
          h_corr = fminf(h_corr, p->h);
          h_corr = fmaxf(h_corr, -p->h * 0.5f);
Pedro Gonnet's avatar
Pedro Gonnet committed
487
        }
488
489

        /* Did we get the right number density? */
490
        if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
491
492
493
494

          /* Ok, correct then */
          p->h += h_corr;

495
          /* Flag for another round of fun */
496
497
          pid[redo] = pid[i];
          redo += 1;
498

499
500
          /* Re-initialise everything */
          hydro_init_part(p);
501

502
          /* Off we go ! */
503
          continue;
504
505
        }

Matthieu Schaller's avatar
Matthieu Schaller committed
506
        /* We now have a particle whose smoothing length has converged */
Matthieu Schaller's avatar
Matthieu Schaller committed
507

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

510
        /* Compute variables required for the force loop */
511
        hydro_prepare_force(p, xp, ti_current, timeBase);
512

Matthieu Schaller's avatar
Matthieu Schaller committed
513
        /* The particle force values are now set.  Do _NOT_
514
           try to read any particle density variables! */
Matthieu Schaller's avatar
Matthieu Schaller committed
515

516
517
        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
518
519
520
      }
    }

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

524
525
526
527
528
529
    /* Re-set the counter for the next loop (potentially). */
    count = redo;
    if (count > 0) {

      /* Climb up the cell hierarchy. */
      for (finger = c; finger != NULL; finger = finger->parent) {
Matthieu Schaller's avatar
Matthieu Schaller committed
530

531
532
        /* Run through this cell's density interactions. */
        for (struct link *l = finger->density; l != NULL; l = l->next) {
Matthieu Schaller's avatar
Matthieu Schaller committed
533

534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
          /* Self-interaction? */
          if (l->t->type == task_type_self)
            runner_doself_subset_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_density(r, finger, parts, pid, count,
                                           l->t->cj);
            else
              runner_dopair_subset_density(r, finger, parts, pid, count,
                                           l->t->ci);

          }

551
552
553
554
555
556
557
          /* 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) {
558
559
560
561
562
563
564
565
566
567
568

            /* 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);
          }
        }
      }
569
    }
570
  }
571

Matthieu Schaller's avatar
Matthieu Schaller committed
572
573
  if (count)
    message("Smoothing length failed to converge on %i particles.", count);
574

Matthieu Schaller's avatar
Matthieu Schaller committed
575
  TIMER_TOC(timer_do_ghost);
576
577
}

578
/**
Pedro Gonnet's avatar
Pedro Gonnet committed
579
 * @brief Mapper function to drift particles and g-particles forward in time.
580
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
581
582
583
 * @param map_data An array of #cells.
 * @param num_elements Chunk size.
 * @param extra_data Pointer to an #engine.
584
 */
585

586
587
588
589
590
591
592
593
594
595
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
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
void runner_do_drift_mapper(void *map_data, int num_elements,
                            void *extra_data) {

  struct engine *e = (struct engine *)extra_data;
  struct cell *cells = (struct cell *)map_data;
  const double timeBase = e->timeBase;
  const double dt = (e->ti_current - e->ti_old) * timeBase;
  const int ti_old = e->ti_old;
  const int ti_current = e->ti_current;

  for (int ind = 0; ind < num_elements; ind++) {
    struct cell *c = &cells[ind];
    if (c == NULL) continue;

    struct part *const parts = c->parts;
    struct xpart *const xparts = c->xparts;
    struct gpart *const gparts = c->gparts;
    float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;

    double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
    double mom[3] = {0.0, 0.0, 0.0};
    double ang_mom[3] = {0.0, 0.0, 0.0};

#ifdef TASK_VERBOSE
    OUT;
#endif

    /* No children? */
    if (!c->split) {

      /* Loop over all the g-particles in the cell */
      const int nr_gparts = c->gcount;
      for (size_t k = 0; k < nr_gparts; k++) {

        /* Get a handle on the gpart. */
        struct gpart *const gp = &gparts[k];

        /* Drift... */
        drift_gpart(gp, dt, timeBase, ti_old, ti_current);

        /* Compute (square of) motion since last cell construction */
        const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
                          gp->x_diff[1] * gp->x_diff[1] +
                          gp->x_diff[2] * gp->x_diff[2];
        dx2_max = fmaxf(dx2_max, dx2);
      }

      /* Loop over all the particles in the cell (more work for these !) */
      const size_t nr_parts = c->count;
      for (size_t k = 0; k < nr_parts; k++) {

        /* Get a handle on the part. */
        struct part *const p = &parts[k];
        struct xpart *const xp = &xparts[k];

        /* Drift... */
        drift_part(p, xp, dt, timeBase, ti_old, ti_current);

        /* Compute (square of) motion since last cell construction */
        const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
                          xp->x_diff[1] * xp->x_diff[1] +
                          xp->x_diff[2] * xp->x_diff[2];
        dx2_max = fmaxf(dx2_max, dx2);

        /* Maximal smoothing length */
        h_max = fmaxf(p->h, h_max);

        /* Now collect quantities for statistics */

        const float half_dt =
            (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
        const double x[3] = {p->x[0], p->x[1], p->x[2]};
        const float v[3] = {xp->v_full[0] + p->a_hydro[0] * half_dt,
                            xp->v_full[1] + p->a_hydro[1] * half_dt,
                            xp->v_full[2] + p->a_hydro[2] * half_dt};
        const float m = p->mass;

        /* Collect mass */
        mass += m;

        /* Collect momentum */
        mom[0] += m * v[0];
        mom[1] += m * v[1];
        mom[2] += m * v[2];

        /* Collect angular momentum */
        ang_mom[0] += m * (x[1] * v[2] - x[2] * v[1]);
        ang_mom[1] += m * (x[2] * v[0] - x[0] * v[2]);
        ang_mom[2] += m * (x[0] * v[1] - x[1] * v[0]);

        /* Collect energies. */
        e_kin += 0.5 * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
        e_pot += 0.;
        e_int += m * hydro_get_internal_energy(p, half_dt);
      }

      /* Now, get the maximal particle motion from its square */
      dx_max = sqrtf(dx2_max);
    }

    /* Otherwise, aggregate data from children. */
    else {

      /* Loop over the progeny and collect their data. */
      for (int k = 0; k < 8; k++)
        if (c->progeny[k] != NULL) {
          struct cell *cp = c->progeny[k];
Pedro Gonnet's avatar
Pedro Gonnet committed
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
          /* Recurse. */
          runner_do_drift_mapper(cp, 1, e);

          dx_max = fmaxf(dx_max, cp->dx_max);
          h_max = fmaxf(h_max, cp->h_max);
          mass += cp->mass;
          e_kin += cp->e_kin;
          e_int += cp->e_int;
          e_pot += cp->e_pot;
          mom[0] += cp->mom[0];
          mom[1] += cp->mom[1];
          mom[2] += cp->mom[2];
          ang_mom[0] += cp->ang_mom[0];
          ang_mom[1] += cp->ang_mom[1];
          ang_mom[2] += cp->ang_mom[2];
        }
    }

    /* Store the values */
    c->h_max = h_max;
    c->dx_max = dx_max;
    c->mass = mass;
    c->e_kin = e_kin;
    c->e_int = e_int;
    c->e_pot = e_pot;
    c->mom[0] = mom[0];
    c->mom[1] = mom[1];
    c->mom[2] = mom[2];
    c->ang_mom[0] = ang_mom[0];
    c->ang_mom[1] = ang_mom[1];
    c->ang_mom[2] = ang_mom[2];
  }
}

728
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
729
730
 * @brief Kick particles in momentum space and collect statistics (fixed
 * time-step case)
731
732
733
 *
 * @param r The runner thread.
 * @param c The cell.
734
 * @param timer Are we timing this ?
735
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
736
void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
737

Matthieu Schaller's avatar
Matthieu Schaller committed
738
  const double global_dt = r->e->dt_max;
739
  const double timeBase = r->e->timeBase;
Matthieu Schaller's avatar
Matthieu Schaller committed
740
  const int count = c->count;
741
742
743
744
  const int gcount = c->gcount;
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
745

746
  int updated = 0, g_updated = 0;
747
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
748
749
750

  TIMER_TIC

Tom Theuns's avatar
Tom Theuns committed
751
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
752
  OUT;
Tom Theuns's avatar
Tom Theuns committed
753
754
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
755
756
757
  /* The new time-step */
  const int new_dti = global_dt / timeBase;

758
759
760
  /* No children? */
  if (!c->split) {

761
    /* Loop over the g-particles and kick everyone. */
762
763
764
765
766
    for (int k = 0; k < gcount; k++) {

      /* Get a handle on the part. */
      struct gpart *const gp = &gparts[k];

767
      /* If the g-particle has no counterpart */
768
      if (gp->id_or_neg_offset > 0) {
769

Matthieu Schaller's avatar
Matthieu Schaller committed
770
771
        /* First, finish the force calculation */
        gravity_end_force(gp);
772

773
774
        /* Kick the g-particle forward */
        kick_gpart(gp, new_dti, timeBase);
775

Matthieu Schaller's avatar
Matthieu Schaller committed
776
777
        /* Number of updated g-particles */
        g_updated++;
778

Matthieu Schaller's avatar
Matthieu Schaller committed
779
780
781
782
783
        /* Minimal time for next end of time-step */
        ti_end_min = min(gp->ti_end, ti_end_min);
        ti_end_max = max(gp->ti_end, ti_end_max);
      }
    }
784

Matthieu Schaller's avatar
Matthieu Schaller committed
785
    /* Now do the hydro ones... */
786

787
    /* Loop over the particles and kick everyone. */
Matthieu Schaller's avatar
Matthieu Schaller committed
788
    for (int k = 0; k < count; k++) {
789

Matthieu Schaller's avatar
Matthieu Schaller committed
790
791
792
      /* Get a handle on the part. */
      struct part *const p = &parts[k];
      struct xpart *const xp = &xparts[k];
793

Matthieu Schaller's avatar
Matthieu Schaller committed
794
795
      /* First, finish the force loop */
      p->h_dt *= p->h * 0.333333333f;
796

Matthieu Schaller's avatar
Matthieu Schaller committed
797
798
799
      /* And do the same of the extra variable */
      hydro_end_force(p);
      if (p->gpart != NULL) gravity_end_force(p->gpart);
800

801
802
      /* Kick the particle forward */
      kick_part(p, xp, new_dti, timeBase);
803

Matthieu Schaller's avatar
Matthieu Schaller committed
804
805
806
      /* Number of updated particles */
      updated++;
      if (p->gpart != NULL) g_updated++;
807

Matthieu Schaller's avatar
Matthieu Schaller committed
808
809
810
      /* Minimal time for next end of time-step */
      ti_end_min = min(p->ti_end, ti_end_min);
      ti_end_max = max(p->ti_end, ti_end_max);
Matthieu Schaller's avatar
Matthieu Schaller committed
811
812
    }
  }
813

Matthieu Schaller's avatar
Matthieu Schaller committed
814
815
  /* Otherwise, aggregate data from children. */
  else {
816

Matthieu Schaller's avatar
Matthieu Schaller committed
817
818
819
820
    /* Loop over the progeny. */
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
        struct cell *const cp = c->progeny[k];
821

Matthieu Schaller's avatar
Matthieu Schaller committed
822
823
        /* Recurse */
        runner_do_kick_fixdt(r, cp, 0);
824

Matthieu Schaller's avatar
Matthieu Schaller committed
825
826
827
828
829
830
831
        /* And aggregate */
        updated += cp->updated;
        g_updated += cp->g_updated;
        ti_end_min = min(cp->ti_end_min, ti_end_min);
        ti_end_max = max(cp->ti_end_max, ti_end_max);
      }
  }
832

Matthieu Schaller's avatar
Matthieu Schaller committed
833
834
835
836
837
  /* Store the values. */
  c->updated = updated;
  c->g_updated = g_updated;
  c->ti_end_min = ti_end_min;
  c->ti_end_max = ti_end_max;
838

Matthieu Schaller's avatar
Matthieu Schaller committed
839
840
  if (timer) TIMER_TOC(timer_kick);
}
841

Matthieu Schaller's avatar
Matthieu Schaller committed
842
843
844
/**
 * @brief Kick particles in momentum space and collect statistics (floating
 * time-step case)
845
846
847
 *
 * @param r The runner thread.
 * @param c The cell.
848
 * @param timer Are we timing this ?
849
 */
850
void runner_do_kick(struct runner *r, struct cell *c, int timer) {
851

852
853
  const struct engine *e = r->e;
  const double timeBase = e->timeBase;
854
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
855
  const int count = c->count;
856
857
858
859
  const int gcount = c->gcount;
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
860

861
  int updated = 0, g_updated = 0;
862
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
863
864
865

  TIMER_TIC

Tom Theuns's avatar
Tom Theuns committed
866
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
867
  OUT;
Tom Theuns's avatar
Tom Theuns committed
868
869
#endif

870
871
872
  /* No children? */
  if (!c->split) {

873
874
875
876
877
878
879
    /* Loop over the g-particles and kick the active ones. */
    for (int k = 0; k < gcount; k++) {

      /* Get a handle on the part. */
      struct gpart *const gp = &gparts[k];

      /* If the g-particle has no counterpart and needs to be kicked */
880
      if (gp->id_or_neg_offset > 0) {
881

Matthieu Schaller's avatar
Matthieu Schaller committed
882
        if (gp->ti_end <= ti_current) {
883

Matthieu Schaller's avatar
Matthieu Schaller committed
884
885
          /* First, finish the force calculation */
          gravity_end_force(gp);
886

887
          /* Compute the next timestep */
888
          const int new_dti = get_gpart_timestep(gp, e);
889

Matthieu Schaller's avatar
Matthieu Schaller committed
890
          /* Now we have a time step, proceed with the kick */
891
          kick_gpart(gp, new_dti, timeBase);
Matthieu Schaller's avatar
Matthieu Schaller committed
892
893
894
895
896
897
898
899

          /* Number of updated g-particles */
          g_updated++;
        }

        /* Minimal time for next end of time-step */
        ti_end_min = min(gp->ti_end, ti_end_min);
        ti_end_max = max(gp->ti_end, ti_end_max);
900
901
902
903
904
      }
    }

    /* Now do the hydro ones... */

905
    /* Loop over the particles and kick the active ones. */
Matthieu Schaller's avatar
Matthieu Schaller committed
906
    for (int k = 0; k < count; k++) {
907
908

      /* Get a handle on the part. */
Matthieu Schaller's avatar
Matthieu Schaller committed
909
910
      struct part *const p = &parts[k];
      struct xpart *const xp = &xparts[k];
911
912

      /* If particle needs to be kicked */
Matthieu Schaller's avatar
Matthieu Schaller committed
913
      if (p->ti_end <= ti_current) {
914
915

        /* First, finish the force loop */
916
        p->h_dt *= p->h * 0.333333333f;
917
918

        /* And do the same of the extra variable */
919
        hydro_end_force(p);
920
        if (p->gpart != NULL) gravity_end_force(p->gpart);
921

Matthieu Schaller's avatar
Matthieu Schaller committed
922
        /* Compute the next timestep (hydro condition) */
923
        const int new_dti = get_part_timestep(p, xp, e);
924

Matthieu Schaller's avatar
Matthieu Schaller committed
925
        /* Now we have a time step, proceed with the kick */
926
        kick_part(p, xp, new_dti, timeBase);
927

Matthieu Schaller's avatar
Matthieu Schaller committed
928
929
        /* Number of updated particles */
        updated++;
930
        if (p->gpart != NULL) g_updated++;
931
932
933
      }

      /* Minimal time for next end of time-step */
934
935
      ti_end_min = min(p->ti_end, ti_end_min);
      ti_end_max = max(p->ti_end, ti_end_max);
936
    }
937
938
  }

939
  /* Otherwise, aggregate data from children. */
940
941
942
  else {

    /* Loop over the progeny. */
943
    for (int k = 0; k < 8; k++)
944
      if (c->progeny[k] != NULL) {
Matthieu Schaller's avatar
Matthieu Schaller committed
945
        struct cell *const cp = c->progeny[k];
946
947

        /* Recurse */
948
        runner_do_kick(r, cp, 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
949

950
        /* And aggregate */
951
        updated += cp->updated;
952
        g_updated += cp->g_updated;
953
954
        ti_end_min = min(cp->ti_end_min, ti_end_min);
        ti_end_max = max(cp->ti_end_max, ti_end_max);
955
956
957
958
      }
  }

  /* Store the values. */
959
  c->updated = updated;
960
  c->g_updated = g_updated;
961
962
  c->ti_end_min = ti_end_min;
  c->ti_end_max = ti_end_max;
963

Peter W. Draper's avatar
Peter W. Draper committed
964
  if (timer) TIMER_TOC(timer_kick);
965
}
966

967
968
969
970
971
972
973
/**
 * @brief Construct the cell properties from the received particles
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
974
void runner_do_recv_cell(struct runner *r, struct cell *c, int timer) {
975
976
977
978
979

  const struct part *const parts = c->parts;
  const struct gpart *const gparts = c->gparts;
  const size_t nr_parts = c->count;
  const size_t nr_gparts = c->gcount;
Matthieu Schaller's avatar
Matthieu Schaller committed
980
  // const int ti_current = r->e->ti_current;
981
982
983
984
985
986
987

  TIMER_TIC;

  int ti_end_min = max_nr_timesteps;
  int ti_end_max = 0;
  float h_max = 0.f;

988
989
990
991
992
993
994
995
996
997
998
999
1000
  /* If this cell is a leaf, collect the particle data. */
  if (!c->split) {

    /* Collect everything... */
    for (size_t k = 0; k < nr_parts; k++) {
      const int ti_end = parts[k].ti_end;
      // if(ti_end < ti_current) error("Received invalid particle !");
      ti_end_min = min(ti_end_min, ti_end);
      ti_end_max = max(ti_end_max, ti_end);
      h_max = fmaxf(h_max, parts[k].h);
    }
    for (size_t k = 0; k < nr_gparts; k++) {
      const int ti_end = gparts[k].ti_end;
For faster browsing, not all history is shown. View entire blame