runner.c 37.3 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. */
150
    struct gpart *restrict gp = &gparts[i];
151
152

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

155
      external_gravity_acceleration(time, potential, constants, gp);
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
 * @brief Calculate change in thermal state of particles induced
 * by radiative cooling and heating.
Stefan Arridge's avatar
Stefan Arridge committed
165
166
167
168
169
170
171
172
 *
 * @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;
173
  struct xpart *restrict xparts = c->xparts;
Stefan Arridge's avatar
Stefan Arridge committed
174
175
  const int count = c->count;
  const int ti_current = r->e->ti_current;
176
  const struct cooling_function_data *cooling_func = r->e->cooling_func;
Stefan Arridge's avatar
Stefan Arridge committed
177
  const struct phys_const *constants = r->e->physical_constants;
178
  const struct UnitSystem *us = r->e->internalUnits;
Stefan Arridge's avatar
Stefan Arridge committed
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
  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
194
  /* Loop over the parts in this cell. */
Stefan Arridge's avatar
Stefan Arridge committed
195
196
197
198
  for (int i = 0; i < count; i++) {

    /* Get a direct pointer on the part. */
    struct part *restrict p = &parts[i];
199
    struct xpart *restrict xp = &xparts[i];
Stefan Arridge's avatar
Stefan Arridge committed
200

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

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

206
      cooling_cool_part(constants, us, cooling_func, p, xp, dt);
Stefan Arridge's avatar
Stefan Arridge committed
207
208
209
210
211
212
    }
  }

  if (timer) TIMER_TOC(timer_do_cooling);
}

Pedro Gonnet's avatar
Pedro Gonnet committed
213
214
215
216
217
218
/**
 * @brief Sort the entries in ascending order using QuickSort.
 *
 * @param sort The entries
 * @param N The number of entries.
 */
219
void runner_do_sort_ascending(struct entry *sort, int N) {
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
271
272
273

  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
274
        }
275
276
277
278
279
280
281
282
283
284
285
286
      } 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
287
    }
288
289
290
  }
}

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

  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
308
309
  float buff[8];
  double px[3];
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333

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

    /* 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
387
        }
388

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

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

429
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
430
  /* Verify the sorting. */
431
432
433
434
435
436
437
438
439
440
  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
441
442
443
444

  if (clock) TIMER_TOC(timer_dosort);
}

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

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

  TIMER_TIC;
461

462
463
464
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

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

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

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

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

480
481
        /* Get ready for a density calculation */
        hydro_init_part(p);
482
      }
483
    }
484
485
486
487
488

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

      /* Get a direct pointer on the part. */
489
      struct gpart *restrict gp = &gparts[i];
490
491
492
493

      if (gp->ti_end <= ti_current) {

        /* Get ready for a density calculation */
494
        gravity_init_gpart(gp);
495
496
      }
    }
497
  }
498

Peter W. Draper's avatar
Peter W. Draper committed
499
  if (timer) TIMER_TOC(timer_init);
500
501
}

502
/**
503
504
505
506
507
508
509
 * @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) {
510

511
#ifdef EXTRA_HYDRO_LOOP
512

513
514
515
516
  struct part *restrict parts = c->parts;
  const int count = c->count;
  const int ti_current = r->e->ti_current;

517
518
519
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
  /* 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);
      }
    }
  }
540
541
542

#else
  error("SWIFT was not compiled with the extra hydro loop activated.");
543
#endif
544
}
545

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

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

568
569
  TIMER_TIC;

570
571
572
  /* Anything to do here? */
  if (c->ti_end_min > ti_current) return;

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

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

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

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

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

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

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

603
        /* Finish the density calculation */
604
        hydro_end_density(p, ti_current);
605

Matthieu Schaller's avatar
Matthieu Schaller committed
606
607
        float h_corr = 0.f;

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

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

          /* Truncate to the range [ -p->h/2 , p->h ]. */
616
617
          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
618
        }
619
620

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

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

626
          /* Flag for another round of fun */
627
628
          pid[redo] = pid[i];
          redo += 1;
629

630
631
          /* Re-initialise everything */
          hydro_init_part(p);
632

633
          /* Off we go ! */
634
          continue;
635
636
        }

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

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

641
        /* Compute variables required for the force loop */
642
        hydro_prepare_force(p, xp, ti_current, timeBase);
643

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

647
648
        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
649
650
651
      }
    }

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

655
656
657
658
659
    /* 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
660
      for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
Matthieu Schaller's avatar
Matthieu Schaller committed
661

662
663
        /* 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
664

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

          }

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

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

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

706
707
708
  /* Be clean */
  free(pid);

Matthieu Schaller's avatar
Matthieu Schaller committed
709
  TIMER_TOC(timer_do_ghost);
710
711
}

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

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

727
728
729
  /* Do we need to drift ? */
  if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return;

730
731
732
  /* Check that we are actually going to move forward. */
  if (ti_current == ti_old) return;

733
734
735
736
  /* 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;
737
738
  double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, e_rad = 0.0;
  double entropy = 0.0, mass = 0.0;
739
740
741
  double mom[3] = {0.0, 0.0, 0.0};
  double ang_mom[3] = {0.0, 0.0, 0.0};

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

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

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

      /* Drift... */
753
      drift_gpart(gp, dt, timeBase, ti_old, ti_current);
754
755
756
757
758

      /* 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];
759
      dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
760
761
762
    }

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

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

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

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

      /* Maximal smoothing length */
780
      h_max = (h_max > p->h) ? h_max : p->h;
781
782
783
784
785
786
787
788
789

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

      const float m = hydro_get_mass(p);
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808

      /* 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.;
809
      e_int += m * hydro_get_internal_energy(p, half_dt);
810
      e_rad += cooling_get_radiated_energy(xp);
811
812

      /* Collect entropy */
813
      entropy += m * hydro_get_entropy(p, half_dt);
814
    }
815
816

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

Matthieu Schaller's avatar
Matthieu Schaller committed
820
821
822
  /* Otherwise, aggregate data from children. */
  else {

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

828
829
        /* Recurse. */
        runner_do_drift(cp, e);
Matthieu Schaller's avatar
Matthieu Schaller committed
830

831
832
        dx_max = max(dx_max, cp->dx_max);
        h_max = max(h_max, cp->h_max);
833
834
835
836
        mass += cp->mass;
        e_kin += cp->e_kin;
        e_int += cp->e_int;
        e_pot += cp->e_pot;
837
        e_rad += cp->e_rad;
838
        entropy += cp->entropy;
839
840
841
842
843
844
        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
845
846
847
848
849
850
      }
  }

  /* Store the values */
  c->h_max = h_max;
  c->dx_max = dx_max;
851
852
853
854
  c->mass = mass;
  c->e_kin = e_kin;
  c->e_int = e_int;
  c->e_pot = e_pot;
855
  c->e_rad = e_rad;
856
  c->entropy = entropy;
857
858
859
860
861
862
  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];
863

864
865
  /* Update the time of the last drift */
  c->ti_old = ti_current;
866
}
867

868
869
870
871
872
873
874
875
876
877
/**
 * @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) {
878

879
880
881
882
883
884
885
  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
886
    if (c != NULL && c->nodeID == e->nodeID) runner_do_drift(c, e);
887
  }
888
}
889

890
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
891
892
 * @brief Kick particles in momentum space and collect statistics (fixed
 * time-step case)
893
894
895
 *
 * @param r The runner thread.
 * @param c The cell.
896
 * @param timer Are we timing this ?
897
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
898
void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
899

Matthieu Schaller's avatar
Matthieu Schaller committed
900
  const double global_dt = r->e->dt_max;
901
  const double timeBase = r->e->timeBase;
Matthieu Schaller's avatar
Matthieu Schaller committed
902
  const int count = c->count;
903
  const int gcount = c->gcount;
Matthieu Schaller's avatar
Matthieu Schaller committed
904
905
906
  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
907
  const double const_G = r->e->physical_constants->const_newton_G;
Matthieu Schaller's avatar
Matthieu Schaller committed
908

909
  int updated = 0, g_updated = 0;
910
  int ti_end_min = max_nr_timesteps, ti_end_max = 0;
911
912
913

  TIMER_TIC

Tom Theuns's avatar
Tom Theuns committed
914
#ifdef TASK_VERBOSE
Matthieu Schaller's avatar
Matthieu Schaller committed
915
  OUT;
Tom Theuns's avatar
Tom Theuns committed
916
917
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
918
919
920
  /* The new time-step */
  const int new_dti = global_dt / timeBase;

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

924
    /* Loop over the g-particles and kick everyone. */
925
926
927
    for (int k = 0; k < gcount; k++) {

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

930
      /* If the g-particle has no counterpart */
931
      if (gp->id_or_neg_offset > 0) {
932

Matthieu Schaller's avatar
Matthieu Schaller committed
933
        /* First, finish the force calculation */
934
        gravity_end_force(gp, const_G);
935

936
937
        /* Kick the g-particle forward */
        kick_gpart(gp, new_dti, timeBase);
938

Matthieu Schaller's avatar
Matthieu Schaller committed
939
940
        /* Number of updated g-particles */
        g_updated++;
941

Matthieu Schaller's avatar
Matthieu Schaller committed
942
943
944
945
946
        /* 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);
      }
    }
947

Matthieu Schaller's avatar
Matthieu Schaller committed
948
    /* Now do the hydro ones... */
949

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

Matthieu Schaller's avatar
Matthieu Schaller committed
953
      /* Get a handle on the part. */
954
955
      struct part *restrict p = &parts[k];
      struct xpart *restrict xp = &xparts[k];
956

Matthieu Schaller's avatar
Matthieu Schaller committed
957
958
      /* First, finish the force loop */
      hydro_end_force(p);
959
      if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
960

961
962
      /* Kick the particle forward */
      kick_part(p, xp, new_dti, timeBase);
963

Matthieu Schaller's avatar
Matthieu Schaller committed
964
965
966
      /* Number of updated particles */
      updated++;
      if (p->gpart != NULL) g_updated++;
967

Matthieu Schaller's avatar
Matthieu Schaller committed
968
969
970
      /* 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
971
972
    }
  }
973

Matthieu Schaller's avatar
Matthieu Schaller committed
974
975
  /* Otherwise, aggregate data from children. */
  else {
976

Matthieu Schaller's avatar
Matthieu Schaller committed
977
978
979
    /* Loop over the progeny. */
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) {
980
        struct cell *restrict cp = c->progeny[k];
981

Matthieu Schaller's avatar
Matthieu Schaller committed
982
983
        /* Recurse */
        runner_do_kick_fixdt(r, cp, 0);
984

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

Matthieu Schaller's avatar
Matthieu Schaller committed
993
994
995
996
997
  /* 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;