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

61
62
63
64
65
66
67
68
69
70
71
72
/**
 * @brief  Entry in a list of sorted indices.
 */
struct entry {

  /*! Distance on the axis */
  float d;

  /*! Particle index */
  int i;
};

73
/* Orientation of the cell pairs */
Matthieu Schaller's avatar
Matthieu Schaller committed
74
75
const double runner_shift[13][3] = {
    {5.773502691896258e-01, 5.773502691896258e-01, 5.773502691896258e-01},
76
    {7.071067811865475e-01, 7.071067811865475e-01, 0.0},
Matthieu Schaller's avatar
Matthieu Schaller committed
77
78
79
80
81
82
83
84
85
86
87
    {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},
Matthieu Schaller's avatar
Matthieu Schaller committed
88
};
89
90

/* Does the axis need flipping ? */
91
92
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
93

94
/* Import the density loop functions. */
95
96
97
#define FUNCTION density
#include "runner_doiact.h"

98
/* Import the gradient loop functions (if required). */
99
100
101
102
103
104
#ifdef EXTRA_HYDRO_LOOP
#undef FUNCTION
#define FUNCTION gradient
#include "runner_doiact.h"
#endif

105
/* Import the force loop functions. */
106
107
108
109
#undef FUNCTION
#define FUNCTION force
#include "runner_doiact.h"

110
/* Import the gravity loop functions. */
111
#include "runner_doiact_fft.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
112
#include "runner_doiact_grav.h"
113

Tom Theuns's avatar
Tom Theuns committed
114
115
116
/**
 * @brief Calculate gravity acceleration from external potential
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
117
118
119
 * @param r runner task
 * @param c cell
 * @param timer 1 if the time is to be recorded.
Tom Theuns's avatar
Tom Theuns committed
120
 */
121
void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
Tom Theuns's avatar
Tom Theuns committed
122

Matthieu Schaller's avatar
Matthieu Schaller committed
123
124
  struct gpart *restrict gparts = c->gparts;
  const int gcount = c->gcount;
125
  const int ti_current = r->e->ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
126
  const struct external_potential *potential = r->e->external_potential;
127
  const struct phys_const *constants = r->e->physical_constants;
128
  const double time = r->e->time;
Matthieu Schaller's avatar
Matthieu Schaller committed
129

130
  TIMER_TIC;
Tom Theuns's avatar
Tom Theuns committed
131

132
133
134
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

Tom Theuns's avatar
Tom Theuns committed
135
136
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
137
    for (int k = 0; k < 8; k++)
138
      if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
Tom Theuns's avatar
Tom Theuns committed
139
140
141
142
    return;
  }

#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
143
  OUT;
Tom Theuns's avatar
Tom Theuns committed
144
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
145

146
  /* Loop over the gparts in this cell. */
Matthieu Schaller's avatar
Matthieu Schaller committed
147
  for (int i = 0; i < gcount; i++) {
148

149
    /* Get a direct pointer on the part. */
Matthieu Schaller's avatar
Matthieu Schaller committed
150
    struct gpart *restrict g = &gparts[i];
151
152
153

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

155
      external_gravity(time, potential, constants, g);
156
    }
157
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
158

159
  if (timer) TIMER_TOC(timer_dograv_external);
Tom Theuns's avatar
Tom Theuns committed
160
161
}

Stefan Arridge's avatar
Stefan Arridge committed
162
163
164
165
166
167
168
169
170
171
172
173
/**
 * @brief Calculate change in entropy from cooling
 *
 * @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) {

  struct part *restrict parts = c->parts;
  const int count = c->count;
  const int ti_current = r->e->ti_current;
Stefan Arridge's avatar
Stefan Arridge committed
174
  const struct cooling_data *cooling = r->e->cooling_data;
Stefan Arridge's avatar
Stefan Arridge committed
175
  const struct phys_const *constants = r->e->physical_constants;
176
  const struct UnitSystem *us = r->e->internalUnits;
Stefan Arridge's avatar
Stefan Arridge committed
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
  const double timeBase = r->e->timeBase;

  TIMER_TIC;

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

#ifdef TASK_VERBOSE
  OUT;
#endif

Stefan Arridge's avatar
Stefan Arridge committed
192
  /* Loop over the parts in this cell. */
Stefan Arridge's avatar
Stefan Arridge committed
193
194
195
196
197
  for (int i = 0; i < count; i++) {

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

Stefan Arridge's avatar
Stefan Arridge committed
198
199
    /* Kick has already updated ti_end, so need to check ti_begin */
    if (p->ti_begin == ti_current) {
200
201
202
203

      const double dt = (p->ti_end - p->ti_begin) * timeBase;

      cooling_cool_part(constants, us, cooling, p, dt);
Stefan Arridge's avatar
Stefan Arridge committed
204
205
206
207
208
209
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
210
211
212
213
214
215
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
216
void runner_do_sort_ascending(struct entry *sort, int N) {
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
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270

  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
271
        }
272
273
274
275
276
277
278
279
280
281
282
283
      } 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
284
    }
285
286
287
  }
}

Pedro Gonnet's avatar
Pedro Gonnet committed
288
289
290
291
292
/**
 * @brief Sort the particles in the given cell along all cardinal directions.
 *
 * @param r The #runner.
 * @param c The #cell.
293
 * @param flags Cell flag.
294
295
 * @param clock Flag indicating whether to record the timing or not, needed
 *      for recursive calls.
Pedro Gonnet's avatar
Pedro Gonnet committed
296
 */
297
void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
298
299
300
301
302
303
304

  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;
Matthieu Schaller's avatar
Matthieu Schaller committed
305
306
  float buff[8];
  double px[3];
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330

  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;
331
      if (missing) runner_do_sort(r, c->progeny[k], missing, 0);
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
    }

    /* 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
384
        }
385

386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
      } /* 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;
Matthieu Schaller's avatar
Matthieu Schaller committed
410
411
412
          sort[j * (count + 1) + k].d = px[0] * runner_shift[j][0] +
                                        px[1] * runner_shift[j][1] +
                                        px[2] * runner_shift[j][2];
413
        }
414
    }
415
416
417
418
419
420

    /* 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;
421
        runner_do_sort_ascending(&sort[j * (count + 1)], count);
422
423
424
425
        c->sorted |= (1 << j);
      }
  }

426
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
427
  /* Verify the sorting. */
428
429
430
431
432
433
434
435
436
437
  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.");
    }
  }
#endif
438
439
440
441

  if (clock) TIMER_TOC(timer_dosort);
}

442
443
444
445
446
/**
 * @brief Initialize the particles before the density calculation
 *
 * @param r The runner thread.
 * @param c The cell.
Matthieu Schaller's avatar
Matthieu Schaller committed
447
 * @param timer 1 if the time is to be recorded.
448
 */
449
void runner_do_init(struct runner *r, struct cell *c, int timer) {
450

Matthieu Schaller's avatar
Matthieu Schaller committed
451
452
  struct part *restrict parts = c->parts;
  struct gpart *restrict gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
453
  const int count = c->count;
454
  const int gcount = c->gcount;
455
  const int ti_current = r->e->ti_current;
456
457

  TIMER_TIC;
458

459
460
461
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

462
463
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
464
    for (int k = 0; k < 8; k++)
465
      if (c->progeny[k] != NULL) runner_do_init(r, c->progeny[k], 0);
466
    return;
467
468
  } else {

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

472
      /* Get a direct pointer on the part. */
473
      struct part *restrict p = &parts[i];
474

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

477
478
        /* Get ready for a density calculation */
        hydro_init_part(p);
479
      }
480
    }
481
482
483
484
485

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

      /* Get a direct pointer on the part. */
486
      struct gpart *restrict gp = &gparts[i];
487
488
489
490

      if (gp->ti_end <= ti_current) {

        /* Get ready for a density calculation */
491
        gravity_init_gpart(gp);
492
493
      }
    }
494
  }
495

Peter W. Draper's avatar
Peter W. Draper committed
496
  if (timer) TIMER_TOC(timer_init);
497
498
}

499
/**
500
501
502
503
504
505
506
 * @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.
 */
void runner_do_extra_ghost(struct runner *r, struct cell *c) {
507

508
#ifdef EXTRA_HYDRO_LOOP
509

510
511
512
513
  struct part *restrict parts = c->parts;
  const int count = c->count;
  const int ti_current = r->e->ti_current;

514
515
516
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k]);
    return;
  } 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];

      if (p->ti_end <= ti_current) {

        /* Get ready for a force calculation */
        hydro_end_gradient(p);
      }
    }
  }
537
538
539

#else
  error("SWIFT was not compiled with the extra hydro loop activated.");
540
#endif
541
}
542

543
/**
544
545
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
546
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
547
 * @param r The runner thread.
548
 * @param c The cell.
549
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
550
void runner_do_ghost(struct runner *r, struct cell *c) {
551

Matthieu Schaller's avatar
Matthieu Schaller committed
552
553
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
554
  int redo, count = c->count;
555
556
  const int ti_current = r->e->ti_current;
  const double timeBase = r->e->timeBase;
557
558
559
560
561
562
563
  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;
564

565
566
  TIMER_TIC;

567
568
569
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

570
571
  /* Recurse? */
  if (c->split) {
Matthieu Schaller's avatar
Matthieu Schaller committed
572
    for (int k = 0; k < 8; k++)
Matthieu Schaller's avatar
Matthieu Schaller committed
573
      if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k]);
574
575
576
577
    return;
  }

  /* Init the IDs that have to be updated. */
578
579
580
  int *pid = NULL;
  if ((pid = malloc(sizeof(int) * count)) == NULL)
    error("Can't allocate memory for pid.");
Matthieu Schaller's avatar
Matthieu Schaller committed
581
  for (int k = 0; k < count; k++) pid[k] = k;
582
583

  /* While there are particles that need to be updated... */
584
  for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
Matthieu Schaller's avatar
Matthieu Schaller committed
585
       num_reruns++) {
586
587
588

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

590
    /* Loop over the parts in this cell. */
Matthieu Schaller's avatar
Matthieu Schaller committed
591
    for (int i = 0; i < count; i++) {
592
593

      /* Get a direct pointer on the part. */
594
595
      struct part *restrict p = &parts[pid[i]];
      struct xpart *restrict xp = &xparts[pid[i]];
596
597

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

600
        /* Finish the density calculation */
601
        hydro_end_density(p, ti_current);
602

Matthieu Schaller's avatar
Matthieu Schaller committed
603
604
        float h_corr = 0.f;

605
        /* If no derivative, double the smoothing length. */
606
        if (p->density.wcount_dh == 0.0f) h_corr = p->h;
607
608
609

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

          /* Truncate to the range [ -p->h/2 , p->h ]. */
613
614
          h_corr = (h_corr < p->h) ? h_corr : p->h;
          h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
Pedro Gonnet's avatar
Pedro Gonnet committed
615
        }
616
617

        /* Did we get the right number density? */
618
        if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
619
620
621
622

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

623
          /* Flag for another round of fun */
624
625
          pid[redo] = pid[i];
          redo += 1;
626

627
628
          /* Re-initialise everything */
          hydro_init_part(p);
629

630
          /* Off we go ! */
631
          continue;
632
633
        }

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

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

638
        /* Compute variables required for the force loop */
639
        hydro_prepare_force(p, xp, ti_current, timeBase);
640

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

644
645
        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
646
647
648
      }
    }

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

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

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

659
660
        /* 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
661

662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
          /* 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);

          }

679
680
681
682
683
684
685
          /* 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) {
686
687
688
689
690
691
692
693
694
695
696

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

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

703
704
705
  /* Be clean */
  free(pid);

Matthieu Schaller's avatar
Matthieu Schaller committed
706
  TIMER_TOC(timer_do_ghost);
707
708
}

709
/**
710
 * @brief Drift particles and g-particles in a cell forward in time
711
712
 *
 * @param c The cell.
713
 * @param e The engine.
714
 */
715
static void runner_do_drift(struct cell *c, struct engine *e) {
716

717
  const double timeBase = e->timeBase;
718
  const int ti_old = c->ti_old;
719
  const int ti_current = e->ti_current;
720
721
722
  struct part *const parts = c->parts;
  struct xpart *const xparts = c->xparts;
  struct gpart *const gparts = c->gparts;
723

724
725
726
  /* Do we need to drift ? */
  if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return;

727
728
729
  /* Check that we are actually going to move forward. */
  if (ti_current == ti_old) return;

730
731
732
733
  /* Drift from the last time the cell was drifted to the current time */
  const double dt = (ti_current - ti_old) * timeBase;

  float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
734
  double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
735
736
737
  double mom[3] = {0.0, 0.0, 0.0};
  double ang_mom[3] = {0.0, 0.0, 0.0};

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

741
    /* Loop over all the g-particles in the cell */
742
    const size_t nr_gparts = c->gcount;
Peter W. Draper's avatar
Peter W. Draper committed
743
744
    for (size_t k = 0; k < nr_gparts; k++) {

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

      /* Drift... */
749
      drift_gpart(gp, dt, timeBase, ti_old, ti_current);
750
751
752
753
754

      /* 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];
755
      dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
756
757
758
    }

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

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

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

769
      /* Compute (square of) motion since last cell construction */
770
771
772
      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];
773
      dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
Matthieu Schaller's avatar
Matthieu Schaller committed
774
775

      /* Maximal smoothing length */
776
      h_max = (h_max > p->h) ? h_max : p->h;
777
778
779
780
781
782
783
784
785

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

      const float m = hydro_get_mass(p);
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804

      /* 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.;
805
      e_int += m * hydro_get_internal_energy(p, half_dt);
806
807

      /* Collect entropy */
808
      entropy += m * hydro_get_entropy(p, half_dt);
809
    }
810
811

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

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

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

823
824
        /* Recurse. */
        runner_do_drift(cp, e);
Matthieu Schaller's avatar
Matthieu Schaller committed
825

826
827
        dx_max = max(dx_max, cp->dx_max);
        h_max = max(h_max, cp->h_max);
828
829
830
831
        mass += cp->mass;
        e_kin += cp->e_kin;
        e_int += cp->e_int;
        e_pot += cp->e_pot;
832
        entropy += cp->entropy;
833
834
835
836
837
838
        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
839
840
841
842
843
844
      }
  }

  /* Store the values */
  c->h_max = h_max;
  c->dx_max = dx_max;
845
846
847
848
  c->mass = mass;
  c->e_kin = e_kin;
  c->e_int = e_int;
  c->e_pot = e_pot;
849
  c->entropy = entropy;
850
851
852
853
854
855
  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];
856

857
858
  /* Update the time of the last drift */
  c->ti_old = ti_current;
859
}
860

861
862
863
864
865
866
867
868
869
870
/**
 * @brief Mapper function to drift particles and g-particles forward in time.
 *
 * @param map_data An array of #cell%s.
 * @param num_elements Chunk size.
 * @param extra_data Pointer to an #engine.
 */

void runner_do_drift_mapper(void *map_data, int num_elements,
                            void *extra_data) {
871

872
873
874
875
876
877
878
  struct engine *e = (struct engine *)extra_data;
  struct cell *cells = (struct cell *)map_data;

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

    /* Only drift local particles. */
Peter W. Draper's avatar
Peter W. Draper committed
879
    if (c != NULL && c->nodeID == e->nodeID) runner_do_drift(c, e);
880
  }
881
}
882

883
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
884
885
 * @brief Kick particles in momentum space and collect statistics (fixed
 * time-step case)
886
887
888
 *
 * @param r The runner thread.
 * @param c The cell.
889
 * @param timer Are we timing this ?
890
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
891
void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
892

Matthieu Schaller's avatar
Matthieu Schaller committed
893
  const double global_dt = r->e->dt_max;
894
  const double timeBase = r->e->timeBase;
Matthieu Schaller's avatar
Matthieu Schaller committed
895
  const int count = c->count;
896
  const int gcount = c->gcount;
Matthieu Schaller's avatar
Matthieu Schaller committed
897
898
899
  struct part *restrict parts = c->parts;
  struct xpart *restrict xparts = c->xparts;
  struct gpart *restrict gparts = c->gparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
900
  const double const_G = r->e->physical_constants->const_newton_G;
Matthieu Schaller's avatar
Matthieu Schaller committed
901

902
  int updated = 0, g_updated = 0;
903
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
904
905
906

  TIMER_TIC

Tom Theuns's avatar
Tom Theuns committed
907
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
908
  OUT;
Tom Theuns's avatar
Tom Theuns committed
909
910
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
911
912
913
  /* The new time-step */
  const int new_dti = global_dt / timeBase;

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

917
    /* Loop over the g-particles and kick everyone. */
918
919
920
    for (int k = 0; k < gcount; k++) {

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

923
      /* If the g-particle has no counterpart */
924
      if (gp->id_or_neg_offset > 0) {
925

Matthieu Schaller's avatar
Matthieu Schaller committed
926
        /* First, finish the force calculation */
927
        gravity_end_force(gp, const_G);
928

929
930
        /* Kick the g-particle forward */
        kick_gpart(gp, new_dti, timeBase);
931

Matthieu Schaller's avatar
Matthieu Schaller committed
932
933
        /* Number of updated g-particles */
        g_updated++;
934

Matthieu Schaller's avatar
Matthieu Schaller committed
935
936
937
938
939
        /* 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);
      }
    }
940

Matthieu Schaller's avatar
Matthieu Schaller committed
941
    /* Now do the hydro ones... */
942

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

Matthieu Schaller's avatar
Matthieu Schaller committed
946
      /* Get a handle on the part. */
947
948
      struct part *restrict p = &parts[k];
      struct xpart *restrict xp = &xparts[k];
949

Matthieu Schaller's avatar
Matthieu Schaller committed
950
951
      /* First, finish the force loop */
      hydro_end_force(p);
952
      if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
953

954
955
      /* Kick the particle forward */
      kick_part(p, xp, new_dti, timeBase);
956

Matthieu Schaller's avatar
Matthieu Schaller committed
957
958
959
      /* Number of updated particles */
      updated++;
      if (p->gpart != NULL) g_updated++;
960

Matthieu Schaller's avatar
Matthieu Schaller committed
961
962
963
      /* 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
964
965
    }
  }
966

Matthieu Schaller's avatar
Matthieu Schaller committed
967
968
  /* Otherwise, aggregate data from children. */
  else {
969

Matthieu Schaller's avatar
Matthieu Schaller committed
970
971
972
    /* Loop over the progeny. */
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
973
        struct cell *restrict cp = c->progeny[k];
974

Matthieu Schaller's avatar
Matthieu Schaller committed
975
976
        /* Recurse */
        runner_do_kick_fixdt(r, cp, 0);
977

Matthieu Schaller's avatar
Matthieu Schaller committed
978
979
980
981
982
983
984
        /* 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);
      }
  }
985

Matthieu Schaller's avatar
Matthieu Schaller committed
986
987
988
989
990
  /* 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;
991

Matthieu Schaller's avatar
Matthieu Schaller committed
992
993
  if (timer) TIMER_TOC(timer_kick);
}
994

Matthieu Schaller's avatar
Matthieu Schaller committed
995
996
997
/**
 * @brief Kick particles in momentum space and collect statistics (floating
 * time-step case)
998
999
1000
 *
 * @param r The runner thread.
 * @param c The cell.
For faster browsing, not all history is shown. View entire blame