runner.c 31.1 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 "const.h"
44
#include "debug.h"
45
#include "drift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
46
#include "engine.h"
47
#include "error.h"
48
49
#include "gravity.h"
#include "hydro.h"
50
#include "hydro_properties.h"
51
#include "kick.h"
52
#include "minmax.h"
53
54
55
56
#include "scheduler.h"
#include "space.h"
#include "task.h"
#include "timers.h"
57
#include "timestep.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
58

59
/* Orientation of the cell pairs */
60
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},
    {0.0, 0.0, 1.0},
74
};
75
76

/* Does the axis need flipping ? */
77
78
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
79

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

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

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

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

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

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

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

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

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

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

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

Pedro Gonnet's avatar
Pedro Gonnet committed
132
133
134
135
136
137
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
138

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

221
void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
222
223
224
225
226
227
228

  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;
229
230
  float buff[8];
  double px[3];
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254

  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;
255
      if (missing) runner_do_sort(r, c->progeny[k], missing, 0);
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
    }

    /* 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
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
      } /* 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;
334
335
336
          sort[j * (count + 1) + k].d = px[0] * runner_shift[j][0] +
                                        px[1] * runner_shift[j][1] +
                                        px[2] * runner_shift[j][2];
337
        }
338
    }
339
340
341
342
343
344

    /* 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;
345
        runner_do_sort_ascending(&sort[j * (count + 1)], count);
346
347
348
349
        c->sorted |= (1 << j);
      }
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
350
351
352
353
354
355
356
357
358
359
360
361
  /* Verify the sorting. */
  /* 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." );
          }
      } */
362
363
364
365

  if (clock) TIMER_TOC(timer_dosort);
}

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

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

  TIMER_TIC;
382

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

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

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

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

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

    /* 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);
      }
    }
415
  }
416

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

420
421
422
/**
 * @brief Intermediate task between density and force
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
423
 * @param r The runner thread.
424
 * @param c The cell.
425
 */
426

Matthieu Schaller's avatar
Matthieu Schaller committed
427
void runner_do_ghost(struct runner *r, struct cell *c) {
428
429

  struct part *p, *parts = c->parts;
430
  struct xpart *xp, *xparts = c->xparts;
431
  struct cell *finger;
Matthieu Schaller's avatar
Matthieu Schaller committed
432
  int redo, count = c->count;
433
  int *pid;
434
  float h_corr;
435
436
  const int ti_current = r->e->ti_current;
  const double timeBase = r->e->timeBase;
437
438
439
440
441
442
443
  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;
444

445
446
  TIMER_TIC;

447
448
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
449
    for (int k = 0; k < 8; k++)
Matthieu Schaller's avatar
Matthieu Schaller committed
450
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k]);
451
452
453
454
455
456
    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
457
  for (int k = 0; k < count; k++) pid[k] = k;
458
459

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

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

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

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

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

476
        /* Finish the density calculation */
477
        hydro_end_density(p, ti_current);
478
479

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

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

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

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

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

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

501
502
          /* Re-initialise everything */
          hydro_init_part(p);
503

504
          /* Off we go ! */
505
          continue;
506
507
        }

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

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

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

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

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

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

526
527
528
529
530
531
    /* 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
532

533
534
        /* 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
535

536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
          /* 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);

          }

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

            /* 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);
          }
        }
      }
571
    }
572
  }
573

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

Matthieu Schaller's avatar
Matthieu Schaller committed
577
  TIMER_TOC(timer_do_ghost);
578
579
}

580
/**
581
 * @brief Drift particles and g-particles forward in time
582
583
584
 *
 * @param r The runner thread.
 * @param c The cell.
585
 * @param timer Are we timing this ?
586
 */
587
void runner_do_drift(struct runner *r, struct cell *c, int timer) {
588

589
590
  const double timeBase = r->e->timeBase;
  const double dt = (r->e->ti_current - r->e->ti_old) * timeBase;
591
592
593
594
595
  const int ti_old = r->e->ti_old;
  const int ti_current = r->e->ti_current;
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;
596
  float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
597

598
599
600
601
  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};

602
  TIMER_TIC
603

Tom Theuns's avatar
Tom Theuns committed
604
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
605
  OUT;
Tom Theuns's avatar
Tom Theuns committed
606
#endif
607

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

611
    /* Loop over all the g-particles in the cell */
Peter W. Draper's avatar
Peter W. Draper committed
612
613
614
    const int nr_gparts = c->gcount;
    for (size_t k = 0; k < nr_gparts; k++) {

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

      /* Drift... */
619
      drift_gpart(gp, dt, timeBase, ti_old, ti_current);
620
621
622
623
624
625

      /* 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);
626
627
628
    }

    /* Loop over all the particles in the cell (more work for these !) */
Peter W. Draper's avatar
Peter W. Draper committed
629
630
    const size_t nr_parts = c->count;
    for (size_t k = 0; k < nr_parts; k++) {
631

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

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

639
      /* Compute (square of) motion since last cell construction */
640
641
642
      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];
643
      dx2_max = fmaxf(dx2_max, dx2);
Matthieu Schaller's avatar
Matthieu Schaller committed
644
645
646

      /* Maximal smoothing length */
      h_max = fmaxf(p->h, h_max);
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

      /* 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.;
674
      e_int += m * hydro_get_internal_energy(p, half_dt);
675
    }
676
677

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

Matthieu Schaller's avatar
Matthieu Schaller committed
681
682
683
684
685
686
  /* Otherwise, aggregate data from children. */
  else {

    /* Loop over the progeny. */
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
687
688

        /* Recurse */
Matthieu Schaller's avatar
Matthieu Schaller committed
689
        struct cell *cp = c->progeny[k];
690
        runner_do_drift(r, cp, 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
691

692
        /* Collect */
693
694
        dx_max = fmaxf(dx_max, cp->dx_max);
        h_max = fmaxf(h_max, cp->h_max);
695
696
697
698
699
700
701
702
703
704
        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];
Matthieu Schaller's avatar
Matthieu Schaller committed
705
706
707
708
709
710
      }
  }

  /* Store the values */
  c->h_max = h_max;
  c->dx_max = dx_max;
711
712
713
714
715
716
717
718
719
720
  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];
721

Peter W. Draper's avatar
Peter W. Draper committed
722
  if (timer) TIMER_TOC(timer_drift);
723
}
724

725
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
 * @brief Kick particles in momentum space and collect statistics (fixed
 * time-step case)
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {

  const double global_dt = r->e->dt_max;
  const double timeBase = r->e->timeBase;
  const int count = c->count;
  const int gcount = c->gcount;
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;

  int updated = 0, g_updated = 0;
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;

  TIMER_TIC

#ifdef TASK_VERBOSE
  OUT;
#endif

  /* The new time-step */
  const int new_dti = global_dt / timeBase;

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

758
    /* Loop over the g-particles and kick everyone. */
Matthieu Schaller's avatar
Matthieu Schaller committed
759
760
761
762
763
    for (int k = 0; k < gcount; k++) {

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

764
      /* If the g-particle has no counterpart */
Matthieu Schaller's avatar
Matthieu Schaller committed
765
766
767
768
769
      if (gp->id < 0) {

        /* First, finish the force calculation */
        gravity_end_force(gp);

770
771
        /* Kick the g-particle forward */
        kick_gpart(gp, new_dti, timeBase);
Matthieu Schaller's avatar
Matthieu Schaller committed
772
773
774
775
776
777
778
779
780
781
782
783

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

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

784
    /* Loop over the particles and kick everyone. */
Matthieu Schaller's avatar
Matthieu Schaller committed
785
786
787
788
789
790
791
792
793
794
795
796
797
    for (int k = 0; k < count; k++) {

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

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

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

798
799
      /* Kick the particle forward */
      kick_part(p, xp, new_dti, timeBase);
Matthieu Schaller's avatar
Matthieu Schaller committed
800
801
802
803
804

      /* Number of updated particles */
      updated++;
      if (p->gpart != NULL) g_updated++;

Matthieu Schaller's avatar
Matthieu Schaller committed
805
806
807
      /* 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
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
    }
  }

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

    /* Loop over the progeny. */
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
        struct cell *const cp = c->progeny[k];

        /* Recurse */
        runner_do_kick_fixdt(r, cp, 0);

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

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

  if (timer) TIMER_TOC(timer_kick);
}

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

849
850
  const struct engine *e = r->e;
  const double timeBase = e->timeBase;
851
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
852
  const int count = c->count;
853
854
855
856
  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
857

858
  int updated = 0, g_updated = 0;
859
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
860
861
862

  TIMER_TIC

Tom Theuns's avatar
Tom Theuns committed
863
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
864
  OUT;
Tom Theuns's avatar
Tom Theuns committed
865
866
#endif

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

870
871
872
873
874
875
876
    /* 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 */
877
      if (gp->id < 0) {
878

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

Matthieu Schaller's avatar
Matthieu Schaller committed
881
882
          /* First, finish the force calculation */
          gravity_end_force(gp);
883

884
          /* Compute the next timestep */
885
          const int new_dti = get_gpart_timestep(gp, e);
886

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

          /* 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);
897
898
899
900
901
      }
    }

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

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

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

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

        /* First, finish the force loop */
913
        p->h_dt *= p->h * 0.333333333f;
914
915

        /* And do the same of the extra variable */
916
        hydro_end_force(p);
917
        if (p->gpart != NULL) gravity_end_force(p->gpart);
918

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

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

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

Matthieu Schaller's avatar
Matthieu Schaller committed
930
931
932
      /* 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);
933
    }
934
935
  }

936
  /* Otherwise, aggregate data from children. */
937
938
939
  else {

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

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

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

  /* Store the values. */
956
  c->updated = updated;
957
  c->g_updated = g_updated;
958
959
  c->ti_end_min = ti_end_min;
  c->ti_end_max = ti_end_max;
960

Peter W. Draper's avatar
Peter W. Draper committed
961
  if (timer) TIMER_TOC(timer_kick);
962
}
963

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

  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
977
  // const int ti_current = r->e->ti_current;
978
979
980
981
982
983
984

  TIMER_TIC;

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

985
  /* Collect everything... */
986
  for (size_t k = 0; k < nr_parts; k++) {
987
    const int ti_end = parts[k].ti_end;
Matthieu Schaller's avatar
Matthieu Schaller committed
988
    // if(ti_end < ti_current) error("Received invalid particle !");
989
990
    ti_end_min = min(ti_end_min, ti_end);
    ti_end_max = max(ti_end_max, ti_end);
991
992
993
    h_max = fmaxf(h_max, parts[k].h);
  }
  for (size_t k = 0; k < nr_gparts; k++) {
994
    const int ti_end = gparts[k].ti_end;
Matthieu Schaller's avatar
Matthieu Schaller committed
995
    // if(ti_end < ti_current) error("Received invalid particle !");
996
997
    ti_end_min = min(ti_end_min, ti_end);
    ti_end_max = max(ti_end_max, ti_end);
998
999
  }

1000
  /* ... and store. */
For faster browsing, not all history is shown. View entire blame