engine.c 173 KB
Newer Older
Pedro Gonnet's avatar
Pedro Gonnet committed
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
Peter W. Draper's avatar
Peter W. Draper committed
5
 *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
6
7
8
 *                    Angus Lepper (angus.lepper@ed.ac.uk)
 *               2016 John A. Regan (john.a.regan@durham.ac.uk)
 *                    Tom Theuns (tom.theuns@durham.ac.uk)
9
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
10
11
12
13
 * 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.
14
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
15
16
17
18
 * 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.
19
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
20
21
 * 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/>.
22
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
23
24
25
26
27
28
29
30
31
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
#include <float.h>
#include <limits.h>
#include <sched.h>
32
#include <stdbool.h>
33
34
35
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
36
#include <sys/stat.h>
37
#include <sys/types.h>
38
#include <unistd.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
39

40
41
/* MPI headers. */
#ifdef WITH_MPI
42

43
#include <mpi.h>
44
45
#endif

Angus Lepper's avatar
Angus Lepper committed
46
47
#ifdef HAVE_LIBNUMA
#include <numa.h>
48
49
#endif

50
/* This object's header. */
Pedro Gonnet's avatar
Pedro Gonnet committed
51
#include "engine.h"
52
53

/* Local headers. */
54
#include "active.h"
55
#include "atomic.h"
56
#include "black_holes_properties.h"
57
#include "cell.h"
58
#include "chemistry.h"
59
#include "clocks.h"
60
#include "cooling.h"
61
#include "cosmology.h"
62
63
#include "cycle.h"
#include "debug.h"
64
#include "entropy_floor.h"
65
#include "equation_of_state.h"
66
#include "error.h"
Alexei Borissov's avatar
Alexei Borissov committed
67
#include "feedback.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
68
#include "fof.h"
69
#include "gravity.h"
70
#include "gravity_cache.h"
71
#include "hydro.h"
lhausamm's avatar
lhausamm committed
72
#include "logger.h"
lhausamm's avatar
lhausamm committed
73
#include "logger_io.h"
74
#include "map.h"
75
#include "memuse.h"
76
#include "minmax.h"
77
#include "mpiuse.h"
78
#include "outputlist.h"
79
#include "parallel_io.h"
80
#include "part.h"
81
#include "partition.h"
James Willis's avatar
James Willis committed
82
#include "profiler.h"
83
#include "proxy.h"
84
#include "restart.h"
85
#include "runner.h"
86
87
#include "serial_io.h"
#include "single_io.h"
88
#include "sort_part.h"
89
#include "star_formation.h"
90
#include "star_formation_logger.h"
Loic Hausammann's avatar
Loic Hausammann committed
91
#include "stars_io.h"
92
#include "statistics.h"
93
#include "timers.h"
94
#include "tools.h"
95
#include "units.h"
96
#include "velociraptor_interface.h"
97
#include "version.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
98

99
100
101
/* Particle cache size. */
#define CACHE_SIZE 512

102
103
104
105
106
const char *engine_policy_names[] = {"none",
                                     "rand",
                                     "steal",
                                     "keep",
                                     "block",
107
                                     "cpu tight",
108
                                     "mpi",
109
                                     "numa affinity",
110
                                     "hydro",
111
112
113
114
115
                                     "self gravity",
                                     "external gravity",
                                     "cosmological integration",
                                     "drift everything",
                                     "reconstruct multi-poles",
116
                                     "temperature",
117
                                     "cooling",
James Willis's avatar
James Willis committed
118
                                     "stars",
119
                                     "structure finding",
120
                                     "star formation",
121
                                     "feedback",
122
                                     "black holes",
123
                                     "fof search",
124
                                     "time-step limiter",
Loic Hausammann's avatar
Loic Hausammann committed
125
126
127
                                     "time-step sync",
                                     "logger"
};
Pedro Gonnet's avatar
Pedro Gonnet committed
128

129
130
131
/** The rank of the engine as a global variable (for messages). */
int engine_rank;

132
/** The current step of the engine as a global variable (for messages). */
133
int engine_current_step;
134

135
136
137
extern int engine_max_parts_per_ghost;
extern int engine_max_sparts_per_ghost;

138
139
140
141
/**
 * @brief Link a density/force task to a cell.
 *
 * @param e The #engine.
142
 * @param l A pointer to the #link, will be modified atomically.
143
144
145
146
 * @param t The #task.
 *
 * @return The new #link pointer.
 */
147
void engine_addlink(struct engine *e, struct link **l, struct task *t) {
148

149
150
151
152
153
154
#ifdef SWIFT_DEBUG_CHECKS
  if (t == NULL) {
    error("Trying to link NULL task.");
  }
#endif

155
  /* Get the next free link. */
156
  const size_t ind = atomic_inc(&e->nr_links);
157
  if (ind >= e->size_links) {
158
159
160
    error(
        "Link table overflow. Increase the value of "
        "`Scheduler:links_per_tasks`.");
161
162
  }
  struct link *res = &e->links[ind];
163

164
  /* Set it atomically. */
165
  res->t = t;
166
  res->next = atomic_swap(l, res);
167
}
168

169
/**
170
 * @brief Repartition the cells amongst the nodes.
171
172
173
 *
 * @param e The #engine.
 */
174
void engine_repartition(struct engine *e) {
175

176
#if defined(WITH_MPI) && (defined(HAVE_PARMETIS) || defined(HAVE_METIS))
177

178
179
  ticks tic = getticks();

180
#ifdef SWIFT_DEBUG_CHECKS
181
182
183
184
  /* Be verbose about this. */
  if (e->nodeID == 0 || e->verbose) message("repartitioning space");
  fflush(stdout);

185
  /* Check that all cells have been drifted to the current time */
186
  space_check_drift_point(e->s, e->ti_current, /*check_multipoles=*/0);
187
188
#endif

189
  /* Clear the repartition flag. */
190
  e->forcerepart = 0;
191

192
  /* Nothing to do if only using a single node. Also avoids METIS
193
   * bug that doesn't handle this case well. */
194
  if (e->nr_nodes == 1) return;
195

196
197
  /* Generate the fixed costs include file. */
  if (e->step > 3 && e->reparttype->trigger <= 1.f) {
Peter W. Draper's avatar
Peter W. Draper committed
198
199
    task_dump_stats("partition_fixed_costs.h", e, /* header = */ 1,
                    /* allranks = */ 1);
200
201
  }

202
  /* Do the repartitioning. */
203
  partition_repartition(e->reparttype, e->nodeID, e->nr_nodes, e->s,
204
                        e->sched.tasks, e->sched.nr_tasks);
205

206
207
  /* Partitioning requires copies of the particles, so we need to reduce the
   * memory in use to the minimum, we can free the sorting indices and the
208
209
   * tasks as these will be regenerated at the next rebuild. Also the foreign
   * particle arrays can go as these will be regenerated in proxy exchange. */
210
211

  /* Sorting indices. */
212
  if (e->s->cells_top != NULL) space_free_cells(e->s);
213

214
215
  /* Task arrays. */
  scheduler_free_tasks(&e->sched);
216

217
218
219
  /* Foreign parts. (no need to nullify the cell pointers as the cells
   * will be regenerated) */
  space_free_foreign_parts(e->s, /*clear_cell_pointers=*/0);
220

221
222
223
  /* Now comes the tricky part: Exchange particles between all nodes.
     This is done in two steps, first allreducing a matrix of
     how many particles go from where to where, then re-allocating
224
     the parts array, and emitting the sends and receives.
225
226
227
     Finally, the space, tasks, and proxies need to be rebuilt. */

  /* Redistribute the particles between the nodes. */
Loic Hausammann's avatar
Loic Hausammann committed
228
  engine_redistribute(e);
229
230
231
232
233
234
235

  /* Make the proxies. */
  engine_makeproxies(e);

  /* Tell the engine it should re-build whenever possible */
  e->forcerebuild = 1;

236
237
238
  /* Flag that a repartition has taken place */
  e->step_props |= engine_step_prop_repartition;

239
240
241
  if (e->verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
242
#else
243
  if (e->reparttype->type != REPART_NONE)
244
    error("SWIFT was not compiled with MPI and METIS or ParMETIS support.");
245
246
247

  /* Clear the repartition flag. */
  e->forcerepart = 0;
248
#endif
249
}
250

251
252
253
254
255
256
257
258
259
/**
 * @brief Decide whether trigger a repartition the cells amongst the nodes.
 *
 * @param e The #engine.
 */
void engine_repartition_trigger(struct engine *e) {

#ifdef WITH_MPI

260
261
  const ticks tic = getticks();

262
263
264
  /* Do nothing if there have not been enough steps since the last repartition
   * as we don't want to repeat this too often or immediately after a
   * repartition step. Also nothing to do when requested. */
265
266
  if (e->step - e->last_repartition >= 2 &&
      e->reparttype->type != REPART_NONE) {
267

268
269
270
271
272
    /* If we have fixed costs available and this is step 2 or we are forcing
     * repartitioning then we do a fixed costs one now. */
    if (e->reparttype->trigger > 1 ||
        (e->step == 2 && e->reparttype->use_fixed_costs)) {

273
      if (e->reparttype->trigger > 1) {
Matthieu Schaller's avatar
Matthieu Schaller committed
274
        if ((e->step % (int)e->reparttype->trigger) == 0) e->forcerepart = 1;
275
276
277
      } else {
        e->forcerepart = 1;
      }
278
      e->reparttype->use_ticks = 0;
279
280
281

    } else {

282
283
284
      /* It is only worth checking the CPU loads when we have processed a
       * significant number of all particles as we require all tasks to have
       * timings. */
285
286
287
288
289
      if ((e->updates > 1 &&
           e->updates >= e->total_nr_parts * e->reparttype->minfrac) ||
          (e->g_updates > 1 &&
           e->g_updates >= e->total_nr_gparts * e->reparttype->minfrac)) {

290
291
292
293
294
295
        /* Should we are use the task timings or fixed costs. */
        if (e->reparttype->use_fixed_costs > 1) {
          e->reparttype->use_ticks = 0;
        } else {
          e->reparttype->use_ticks = 1;
        }
296

297
        /* Get CPU time used since the last call to this function. */
Peter W. Draper's avatar
Peter W. Draper committed
298
299
        double elapsed_cputime =
            clocks_get_cputime_used() - e->cputime_last_step;
300
301
302

        /* Gather the elapsed CPU times from all ranks for the last step. */
        double elapsed_cputimes[e->nr_nodes];
Peter W. Draper's avatar
Peter W. Draper committed
303
304
        MPI_Gather(&elapsed_cputime, 1, MPI_DOUBLE, elapsed_cputimes, 1,
                   MPI_DOUBLE, 0, MPI_COMM_WORLD);
305
306
        if (e->nodeID == 0) {

307
          /* Get the range and mean of cputimes. */
308
309
          double mintime = elapsed_cputimes[0];
          double maxtime = elapsed_cputimes[0];
310
          double sum = elapsed_cputimes[0];
311
312
313
          for (int k = 1; k < e->nr_nodes; k++) {
            if (elapsed_cputimes[k] > maxtime) maxtime = elapsed_cputimes[k];
            if (elapsed_cputimes[k] < mintime) mintime = elapsed_cputimes[k];
314
            sum += elapsed_cputimes[k];
315
          }
316
          double mean = sum / (double)e->nr_nodes;
317
318

          /* Are we out of balance? */
319
320
          double abs_trigger = fabs(e->reparttype->trigger);
          if (((maxtime - mintime) / mean) > abs_trigger) {
321
            if (e->verbose)
Peter W. Draper's avatar
Peter W. Draper committed
322
323
              message("trigger fraction %.3f > %.3f will repartition",
                      (maxtime - mintime) / mean, abs_trigger);
324
            e->forcerepart = 1;
325
          } else {
Peter W. Draper's avatar
Peter W. Draper committed
326
            if (e->verbose)
Peter W. Draper's avatar
Peter W. Draper committed
327
328
              message("trigger fraction %.3f =< %.3f will not repartition",
                      (maxtime - mintime) / mean, abs_trigger);
329
330
331
          }
        }
      }
332
333
334

      /* All nodes do this together. */
      MPI_Bcast(&e->forcerepart, 1, MPI_INT, 0, MPI_COMM_WORLD);
335
336
337
    }

    /* Remember we did this. */
Peter W. Draper's avatar
Peter W. Draper committed
338
    if (e->forcerepart) e->last_repartition = e->step;
339
340
  }

341
342
343
  /* We always reset CPU time for next check, unless it will not be used. */
  if (e->reparttype->type != REPART_NONE)
    e->cputime_last_step = clocks_get_cputime_used();
344
345
346
347

  if (e->verbose)
    message("took %.3f %s", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
348
349
350
#endif
}

351
/**
352
353
354
355
 * @brief Exchange cell structures with other nodes.
 *
 * @param e The #engine.
 */
356
void engine_exchange_cells(struct engine *e) {
357

358
#ifdef WITH_MPI
359

360
  const int with_gravity = e->policy & engine_policy_self_gravity;
361
  const ticks tic = getticks();
362

363
  /* Exchange the cell structure with neighbouring ranks. */
364
  proxy_cells_exchange(e->proxies, e->nr_proxies, e->s, with_gravity);
365

366
367
368
  if (e->verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
369
370
371
372
373
374
375

#else
  error("SWIFT was not compiled with MPI support.");
#endif
}

/**
376
 * @brief Exchange straying particles with other nodes.
377
378
 *
 * @param e The #engine.
Pedro Gonnet's avatar
Pedro Gonnet committed
379
 * @param offset_parts The index in the parts array as of which the foreign
380
 *        parts reside (i.e. the current number of local #part).
381
 * @param ind_part The foreign #cell ID of each part.
382
 * @param Npart The number of stray parts, contains the number of parts received
383
 *        on return.
Pedro Gonnet's avatar
Pedro Gonnet committed
384
 * @param offset_gparts The index in the gparts array as of which the foreign
385
 *        parts reside (i.e. the current number of local #gpart).
386
 * @param ind_gpart The foreign #cell ID of each gpart.
Pedro Gonnet's avatar
Pedro Gonnet committed
387
 * @param Ngpart The number of stray gparts, contains the number of gparts
Pedro Gonnet's avatar
Pedro Gonnet committed
388
 *        received on return.
389
 * @param offset_sparts The index in the sparts array as of which the foreign
390
 *        parts reside (i.e. the current number of local #spart).
391
392
393
 * @param ind_spart The foreign #cell ID of each spart.
 * @param Nspart The number of stray sparts, contains the number of sparts
 *        received on return.
394
395
396
397
398
 * @param offset_bparts The index in the bparts array as of which the foreign
 *        parts reside (i.e. the current number of local #bpart).
 * @param ind_bpart The foreign #cell ID of each bpart.
 * @param Nbpart The number of stray bparts, contains the number of bparts
 *        received on return.
399
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
400
401
 * Note that this function does not mess-up the linkage between parts and
 * gparts, i.e. the received particles have correct linkeage.
402
 */
403
404
405
406
void engine_exchange_strays(struct engine *e, const size_t offset_parts,
                            const int *ind_part, size_t *Npart,
                            const size_t offset_gparts, const int *ind_gpart,
                            size_t *Ngpart, const size_t offset_sparts,
407
408
409
                            const int *ind_spart, size_t *Nspart,
                            const size_t offset_bparts, const int *ind_bpart,
                            size_t *Nbpart) {
410
411

#ifdef WITH_MPI
412
  struct space *s = e->s;
413
  ticks tic = getticks();
414
415

  /* Re-set the proxies. */
Pedro Gonnet's avatar
Pedro Gonnet committed
416
417
418
  for (int k = 0; k < e->nr_proxies; k++) {
    e->proxies[k].nr_parts_out = 0;
    e->proxies[k].nr_gparts_out = 0;
419
    e->proxies[k].nr_sparts_out = 0;
420
    e->proxies[k].nr_bparts_out = 0;
Pedro Gonnet's avatar
Pedro Gonnet committed
421
  }
422

423
  /* Put the parts into the corresponding proxies. */
424
  for (size_t k = 0; k < *Npart; k++) {
425
426
427
428

    /* Ignore the particles we want to get rid of (inhibited, ...). */
    if (ind_part[k] == -1) continue;

429
    /* Get the target node and proxy ID. */
430
    const int node_id = e->s->cells_top[ind_part[k]].nodeID;
431
432
    if (node_id < 0 || node_id >= e->nr_nodes)
      error("Bad node ID %i.", node_id);
433
    const int pid = e->proxy_ind[node_id];
434
    if (pid < 0) {
435
436
      error(
          "Do not have a proxy for the requested nodeID %i for part with "
437
          "id=%lld, x=[%e,%e,%e].",
Pedro Gonnet's avatar
Pedro Gonnet committed
438
439
440
          node_id, s->parts[offset_parts + k].id,
          s->parts[offset_parts + k].x[0], s->parts[offset_parts + k].x[1],
          s->parts[offset_parts + k].x[2]);
441
    }
Pedro Gonnet's avatar
Pedro Gonnet committed
442

443
444
    /* Re-link the associated gpart with the buffer offset of the part. */
    if (s->parts[offset_parts + k].gpart != NULL) {
445
446
      s->parts[offset_parts + k].gpart->id_or_neg_offset =
          -e->proxies[pid].nr_parts_out;
447
    }
Pedro Gonnet's avatar
Pedro Gonnet committed
448

449
450
451
452
453
#ifdef SWIFT_DEBUG_CHECKS
    if (s->parts[offset_parts + k].time_bin == time_bin_inhibited)
      error("Attempting to exchange an inhibited particle");
#endif

454
    /* Load the part and xpart into the proxy. */
455
456
    proxy_parts_load(&e->proxies[pid], &s->parts[offset_parts + k],
                     &s->xparts[offset_parts + k], 1);
457
458

#ifdef WITH_LOGGER
Loic Hausammann's avatar
Loic Hausammann committed
459
460
461
462
463
464
465
466
467
468
469
    if (e->policy & engine_policy_logger) {
      const int logger_flag = logger_generate_flag(
        logger_flag_mpi_exit, node_id);

      /* Log the particle when leaving a rank. */
      logger_log_part(e->logger, &s->parts[offset_parts + k],
                      logger_masks_all_part |
                      logger_mask_data[logger_special_flags].mask,
                      &s->xparts[offset_parts + k].logger_data.last_offset,
                      logger_flag);
    }
470
#endif
471
  }
472

473
  /* Put the sparts into the corresponding proxies. */
474
  for (size_t k = 0; k < *Nspart; k++) {
475
476
477
478
479

    /* Ignore the particles we want to get rid of (inhibited, ...). */
    if (ind_spart[k] == -1) continue;

    /* Get the target node and proxy ID. */
480
481
482
483
    const int node_id = e->s->cells_top[ind_spart[k]].nodeID;
    if (node_id < 0 || node_id >= e->nr_nodes)
      error("Bad node ID %i.", node_id);
    const int pid = e->proxy_ind[node_id];
484
    if (pid < 0) {
485
486
487
488
489
490
      error(
          "Do not have a proxy for the requested nodeID %i for part with "
          "id=%lld, x=[%e,%e,%e].",
          node_id, s->sparts[offset_sparts + k].id,
          s->sparts[offset_sparts + k].x[0], s->sparts[offset_sparts + k].x[1],
          s->sparts[offset_sparts + k].x[2]);
491
    }
492

493
    /* Re-link the associated gpart with the buffer offset of the spart. */
494
495
496
497
498
    if (s->sparts[offset_sparts + k].gpart != NULL) {
      s->sparts[offset_sparts + k].gpart->id_or_neg_offset =
          -e->proxies[pid].nr_sparts_out;
    }

499
500
501
502
503
#ifdef SWIFT_DEBUG_CHECKS
    if (s->sparts[offset_sparts + k].time_bin == time_bin_inhibited)
      error("Attempting to exchange an inhibited particle");
#endif

504
505
    /* Load the spart into the proxy */
    proxy_sparts_load(&e->proxies[pid], &s->sparts[offset_sparts + k], 1);
506
507

#ifdef WITH_LOGGER
Loic Hausammann's avatar
Loic Hausammann committed
508
509
510
    if (e->policy & engine_policy_logger) {
      const int logger_flag = logger_generate_flag(
        logger_flag_mpi_exit, node_id);
Loic Hausammann's avatar
Loic Hausammann committed
511
512

      /* Log the particle when leaving a rank. */
Loic Hausammann's avatar
Loic Hausammann committed
513
514
515
516
517
518
      logger_log_spart(e->logger, &s->sparts[offset_sparts + k],
                       logger_masks_all_spart |
                       logger_mask_data[logger_special_flags].mask,
                       &s->sparts[offset_parts + k].logger_data.last_offset,
                       logger_flag);
    }
519
#endif
520
  }
521

522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
  /* Put the bparts into the corresponding proxies. */
  for (size_t k = 0; k < *Nbpart; k++) {

    /* Ignore the particles we want to get rid of (inhibited, ...). */
    if (ind_bpart[k] == -1) continue;

    /* Get the target node and proxy ID. */
    const int node_id = e->s->cells_top[ind_bpart[k]].nodeID;
    if (node_id < 0 || node_id >= e->nr_nodes)
      error("Bad node ID %i.", node_id);
    const int pid = e->proxy_ind[node_id];
    if (pid < 0) {
      error(
          "Do not have a proxy for the requested nodeID %i for part with "
          "id=%lld, x=[%e,%e,%e].",
          node_id, s->bparts[offset_bparts + k].id,
          s->bparts[offset_bparts + k].x[0], s->bparts[offset_bparts + k].x[1],
          s->bparts[offset_bparts + k].x[2]);
    }

    /* Re-link the associated gpart with the buffer offset of the bpart. */
    if (s->bparts[offset_bparts + k].gpart != NULL) {
      s->bparts[offset_bparts + k].gpart->id_or_neg_offset =
          -e->proxies[pid].nr_bparts_out;
    }

#ifdef SWIFT_DEBUG_CHECKS
    if (s->bparts[offset_bparts + k].time_bin == time_bin_inhibited)
      error("Attempting to exchange an inhibited particle");
#endif

    /* Load the bpart into the proxy */
    proxy_bparts_load(&e->proxies[pid], &s->bparts[offset_bparts + k], 1);
555
556

#ifdef WITH_LOGGER
Loic Hausammann's avatar
Loic Hausammann committed
557
558
559
    if (e->policy & engine_policy_logger) {
      error("TODO");
    }
560
#endif
561
  }
562

563
564
  /* Put the gparts into the corresponding proxies. */
  for (size_t k = 0; k < *Ngpart; k++) {
565
566
567
568
569

    /* Ignore the particles we want to get rid of (inhibited, ...). */
    if (ind_gpart[k] == -1) continue;

    /* Get the target node and proxy ID. */
570
571
572
573
    const int node_id = e->s->cells_top[ind_gpart[k]].nodeID;
    if (node_id < 0 || node_id >= e->nr_nodes)
      error("Bad node ID %i.", node_id);
    const int pid = e->proxy_ind[node_id];
574
    if (pid < 0) {
575
576
577
578
579
580
      error(
          "Do not have a proxy for the requested nodeID %i for part with "
          "id=%lli, x=[%e,%e,%e].",
          node_id, s->gparts[offset_gparts + k].id_or_neg_offset,
          s->gparts[offset_gparts + k].x[0], s->gparts[offset_gparts + k].x[1],
          s->gparts[offset_gparts + k].x[2]);
581
582
583
584
585
586
    }

#ifdef SWIFT_DEBUG_CHECKS
    if (s->gparts[offset_gparts + k].time_bin == time_bin_inhibited)
      error("Attempting to exchange an inhibited particle");
#endif
587
588
589

    /* Load the gpart into the proxy */
    proxy_gparts_load(&e->proxies[pid], &s->gparts[offset_gparts + k], 1);
590
591
592

#ifdef WITH_LOGGER
    /* Write only the dark matter particles */
Loic Hausammann's avatar
Loic Hausammann committed
593
594
    if ((e->policy & engine_policy_logger) &&
        s->gparts[offset_gparts + k].type == swift_type_dark_matter) {
Loic Hausammann's avatar
Loic Hausammann committed
595
596

      const int logger_flag = logger_generate_flag(
Loic Hausammann's avatar
Loic Hausammann committed
597
         logger_flag_mpi_exit, node_id);
598
599

      /* Log the particle when leaving a rank. */
Loic Hausammann's avatar
Loic Hausammann committed
600
      logger_log_gpart(e->logger, &s->gparts[offset_gparts + k],
Loic Hausammann's avatar
Loic Hausammann committed
601
                       logger_masks_all_gpart |
602
                       logger_mask_data[logger_special_flags].mask,
Loic Hausammann's avatar
Loic Hausammann committed
603
604
                       &s->sparts[offset_parts + k].logger_data.last_offset,
                       logger_flag);
605
606
    }
#endif
607
608
  }

609
  /* Launch the proxies. */
610
611
  MPI_Request reqs_in[5 * engine_maxproxies];
  MPI_Request reqs_out[5 * engine_maxproxies];
612
  for (int k = 0; k < e->nr_proxies; k++) {
613
    proxy_parts_exchange_first(&e->proxies[k]);
614
615
616
617
618
    reqs_in[k] = e->proxies[k].req_parts_count_in;
    reqs_out[k] = e->proxies[k].req_parts_count_out;
  }

  /* Wait for each count to come in and start the recv. */
619
  for (int k = 0; k < e->nr_proxies; k++) {
620
    int pid = MPI_UNDEFINED;
621
622
    if (MPI_Waitany(e->nr_proxies, reqs_in, &pid, MPI_STATUS_IGNORE) !=
            MPI_SUCCESS ||
623
624
625
        pid == MPI_UNDEFINED)
      error("MPI_Waitany failed.");
    // message( "request from proxy %i has arrived." , pid );
626
    proxy_parts_exchange_second(&e->proxies[pid]);
627
  }
628

629
  /* Wait for all the sends to have finished too. */
630
631
632
  if (MPI_Waitall(e->nr_proxies, reqs_out, MPI_STATUSES_IGNORE) != MPI_SUCCESS)
    error("MPI_Waitall on sends failed.");

633
  /* Count the total number of incoming particles and make sure we have
634
     enough space to accommodate them. */
635
636
  int count_parts_in = 0;
  int count_gparts_in = 0;
637
  int count_sparts_in = 0;
638
  int count_bparts_in = 0;
639
640
641
  for (int k = 0; k < e->nr_proxies; k++) {
    count_parts_in += e->proxies[k].nr_parts_in;
    count_gparts_in += e->proxies[k].nr_gparts_in;
642
    count_sparts_in += e->proxies[k].nr_sparts_in;
643
    count_bparts_in += e->proxies[k].nr_bparts_in;
644
645
  }
  if (e->verbose) {
646
647
648
649
650
    message(
        "sent out %zu/%zu/%zu/%zu parts/gparts/sparts/bparts, got %i/%i/%i/%i "
        "back.",
        *Npart, *Ngpart, *Nspart, *Nbpart, count_parts_in, count_gparts_in,
        count_sparts_in, count_bparts_in);
651
  }
652
653

  /* Reallocate the particle arrays if necessary */
654
  if (offset_parts + count_parts_in > s->size_parts) {
655
    s->size_parts = (offset_parts + count_parts_in) * engine_parts_size_grow;
656
657
    struct part *parts_new = NULL;
    struct xpart *xparts_new = NULL;
658
    if (swift_memalign("parts", (void **)&parts_new, part_align,
659
                       sizeof(struct part) * s->size_parts) != 0 ||
660
        swift_memalign("xparts", (void **)&xparts_new, xpart_align,
661
662
                       sizeof(struct xpart) * s->size_parts) != 0)
      error("Failed to allocate new part data.");
663
664
    memcpy(parts_new, s->parts, sizeof(struct part) * offset_parts);
    memcpy(xparts_new, s->xparts, sizeof(struct xpart) * offset_parts);
665
666
    swift_free("parts", s->parts);
    swift_free("xparts", s->xparts);
667
668
    s->parts = parts_new;
    s->xparts = xparts_new;
669
670

    /* Reset the links */
Pedro Gonnet's avatar
Pedro Gonnet committed
671
672
    for (size_t k = 0; k < offset_parts; k++) {
      if (s->parts[k].gpart != NULL) {
673
        s->parts[k].gpart->id_or_neg_offset = -k;
Pedro Gonnet's avatar
Pedro Gonnet committed
674
675
      }
    }
676
  }
677

678
679
680
  if (offset_sparts + count_sparts_in > s->size_sparts) {
    s->size_sparts = (offset_sparts + count_sparts_in) * engine_parts_size_grow;
    struct spart *sparts_new = NULL;
681
    if (swift_memalign("sparts", (void **)&sparts_new, spart_align,
682
683
684
                       sizeof(struct spart) * s->size_sparts) != 0)
      error("Failed to allocate new spart data.");
    memcpy(sparts_new, s->sparts, sizeof(struct spart) * offset_sparts);
685
    swift_free("sparts", s->sparts);
686
    s->sparts = sparts_new;
687
688

    /* Reset the links */
689
    for (size_t k = 0; k < offset_sparts; k++) {
690
691
692
693
694
      if (s->sparts[k].gpart != NULL) {
        s->sparts[k].gpart->id_or_neg_offset = -k;
      }
    }
  }
695

696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
  if (offset_bparts + count_bparts_in > s->size_bparts) {
    s->size_bparts = (offset_bparts + count_bparts_in) * engine_parts_size_grow;
    struct bpart *bparts_new = NULL;
    if (swift_memalign("bparts", (void **)&bparts_new, bpart_align,
                       sizeof(struct bpart) * s->size_bparts) != 0)
      error("Failed to allocate new bpart data.");
    memcpy(bparts_new, s->bparts, sizeof(struct bpart) * offset_bparts);
    swift_free("bparts", s->bparts);
    s->bparts = bparts_new;

    /* Reset the links */
    for (size_t k = 0; k < offset_bparts; k++) {
      if (s->bparts[k].gpart != NULL) {
        s->bparts[k].gpart->id_or_neg_offset = -k;
      }
    }
  }
713

714
  if (offset_gparts + count_gparts_in > s->size_gparts) {
715
    s->size_gparts = (offset_gparts + count_gparts_in) * engine_parts_size_grow;
716
    struct gpart *gparts_new = NULL;
717
    if (swift_memalign("gparts", (void **)&gparts_new, gpart_align,
718
719
720
                       sizeof(struct gpart) * s->size_gparts) != 0)
      error("Failed to allocate new gpart data.");
    memcpy(gparts_new, s->gparts, sizeof(struct gpart) * offset_gparts);
721
    swift_free("gparts", s->gparts);
722
    s->gparts = gparts_new;
723

724
    /* Reset the links */
Pedro Gonnet's avatar
Pedro Gonnet committed
725
    for (size_t k = 0; k < offset_gparts; k++) {
726
      if (s->gparts[k].type == swift_type_gas) {
727
        s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
728
      } else if (s->gparts[k].type == swift_type_stars) {
729
        s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
730
731
      } else if (s->gparts[k].type == swift_type_black_hole) {
        s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
Pedro Gonnet's avatar
Pedro Gonnet committed
732
733
      }
    }
734
  }
735
736

  /* Collect the requests for the particle data from the proxies. */
737
738
  int nr_in = 0, nr_out = 0;
  for (int k = 0; k < e->nr_proxies; k++) {
739
    if (e->proxies[k].nr_parts_in > 0) {
740
741
      reqs_in[5 * k] = e->proxies[k].req_parts_in;
      reqs_in[5 * k + 1] = e->proxies[k].req_xparts_in;
742
743
      nr_in += 2;
    } else {
744
      reqs_in[5 * k] = reqs_in[5 * k + 1] = MPI_REQUEST_NULL;
745
746
    }
    if (e->proxies[k].nr_gparts_in > 0) {
747
      reqs_in[5 * k + 2] = e->proxies[k].req_gparts_in;
748
      nr_in += 1;
749
    } else {
750
      reqs_in[5 * k + 2] = MPI_REQUEST_NULL;
751
    }
752
    if (e->proxies[k].nr_sparts_in > 0) {
753
      reqs_in[5 * k + 3] = e->proxies[k].req_sparts_in;
754
755
      nr_in += 1;
    } else {
756
757
758
759
760
761
762
      reqs_in[5 * k + 3] = MPI_REQUEST_NULL;
    }
    if (e->proxies[k].nr_bparts_in > 0) {
      reqs_in[5 * k + 4] = e->proxies[k].req_bparts_in;
      nr_in += 1;
    } else {
      reqs_in[5 * k + 4] = MPI_REQUEST_NULL;
763
764
    }

765
    if (e->proxies[k].nr_parts_out > 0) {
766
767
      reqs_out[5 * k] = e->proxies[k].req_parts_out;
      reqs_out[5 * k + 1] = e->proxies[k].req_xparts_out;
768
769
      nr_out += 2;
    } else {
770
      reqs_out[5 * k] = reqs_out[5 * k + 1] = MPI_REQUEST_NULL;
771
772
    }
    if (e->proxies[k].nr_gparts_out > 0) {
773
      reqs_out[5 * k + 2] = e->proxies[k].req_gparts_out;
774
775
      nr_out += 1;
    } else {
776
      reqs_out[5 * k + 2] = MPI_REQUEST_NULL;
777
778
    }
    if (e->proxies[k].nr_sparts_out > 0) {
779
780
781
782
783
784
785
      reqs_out[5 * k + 3] = e->proxies[k].req_sparts_out;
      nr_out += 1;
    } else {
      reqs_out[5 * k + 3] = MPI_REQUEST_NULL;
    }
    if (e->proxies[k].nr_bparts_out > 0) {
      reqs_out[5 * k + 4] = e->proxies[k].req_bparts_out;
786
      nr_out += 1;
787
    } else {
788
      reqs_out[5 * k + 4] = MPI_REQUEST_NULL;
789
    }
790
791
792
793
  }

  /* Wait for each part array to come in and collect the new
     parts from the proxies. */
794
  int count_parts = 0, count_gparts = 0, count_sparts = 0, count_bparts = 0;
795
  for (int k = 0; k < nr_in; k++) {
796
    int err, pid;
797
    if ((err = MPI_Waitany(5 * e->nr_proxies, reqs_in, &pid,
798
                           MPI_STATUS_IGNORE)) != MPI_SUCCESS) {
799
800
801
802
803
804
      char buff[MPI_MAX_ERROR_STRING];
      int res;
      MPI_Error_string(err, buff, &res);
      error("MPI_Waitany failed (%s).", buff);
    }
    if (pid == MPI_UNDEFINED) break;
805
806
    // message( "request from proxy %i has arrived." , pid / 5 );
    pid = 5 * (pid / 5);
Pedro Gonnet's avatar
Pedro Gonnet committed
807

808
    /* If all the requests for a given proxy have arrived... */
809
810
    if (reqs_in[pid + 0] == MPI_REQUEST_NULL &&
        reqs_in[pid + 1] == MPI_REQUEST_NULL &&
811
        reqs_in[pid + 2] == MPI_REQUEST_NULL &&
812
        reqs_in[pid + 3] == MPI_REQUEST_NULL &&
813
        reqs_in[pid + 4] == MPI_REQUEST_NULL) {
814
      /* Copy the particle data to the part/xpart/gpart arrays. */
815
      struct proxy *prox = &e->proxies[pid / 5];
816
817
818
819
820
821
      memcpy(&s->parts[offset_parts + count_parts], prox->parts_in,
             sizeof(struct part) * prox->nr_parts_in);
      memcpy(&s->xparts[offset_parts + count_parts], prox->xparts_in,
             sizeof(struct xpart) * prox->nr_parts_in);
      memcpy(&s->gparts[offset_gparts + count_gparts], prox->gparts_in,
             sizeof(struct gpart) * prox->nr_gparts_in);
822
823
      memcpy(&s->sparts[offset_sparts + count_sparts], prox->sparts_in,
             sizeof(struct spart) * prox->nr_sparts_in);
824
825
      memcpy(&s->bparts[offset_bparts + count_bparts], prox->bparts_in,
             sizeof(struct bpart) * prox->nr_bparts_in);
826
827

#ifdef WITH_LOGGER
Loic Hausammann's avatar
Loic Hausammann committed
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
      if (e->policy & engine_policy_logger) {
        const int flag = logger_generate_flag(logger_flag_mpi_enter,
                                              prox->nodeID);

        struct part *parts = &s->parts[offset_parts + count_parts];
        struct xpart *xparts = &s->xparts[offset_parts + count_parts];
        struct spart *sparts = &s->sparts[offset_sparts + count_sparts];
        struct gpart *gparts = &s->gparts[offset_gparts + count_gparts];

        /* Log the gas particles */
        const unsigned int mask_hydro =
          logger_masks_all_part |
          logger_mask_data[logger_special_flags].mask;

        for(int i = 0; i < prox->nr_parts_in; i++) {
          logger_log_part(e->logger, &parts[i], mask_hydro,
                          &xparts[i].logger_data.last_offset,
                          flag);
          /* Reset the counter */
          xparts[i].logger_data.steps_since_last_output = 0;
        }

        /* Log the stellar particles */
        const unsigned int mask_stars = logger_masks_all_spart |
          logger_mask_data[logger_special_flags].mask;
        for(int i = 0; i < prox->nr_sparts_in; i++) {
          logger_log_spart(e->logger, &sparts[i], mask_stars,
                           &sparts[i].logger_data.last_offset,
                           flag);
          sparts[i].logger_data.steps_since_last_output = 0;
        }

        /* Log the gparts */
        const unsigned int mask_grav =
          logger_masks_all_gpart |
          logger_mask_data[logger_special_flags].mask;
        for(int i = 0; i < prox->nr_gparts_in; i++) {
          /* Log only the dark matter */
          if (gparts[i].type != swift_type_dark_matter) continue;

          logger_log_gpart(e->logger, &gparts[i], mask_grav,
                           &gparts[i].logger_data.last_offset,
                           flag);
          gparts[i].logger_data.steps_since_last_output = 0;
        }

        /* Log the bparts */
        if (prox->nr_bparts_in > 0) {
          error("TODO");
        }
      }
879
#endif
880
881
      /* for (int k = offset; k < offset + count; k++)
         message(
882
            "received particle %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i.",
883
            s->parts[k].id, s->parts[k].x[0], s->parts[k].x[1],
884
            s->parts[k].x[2], s->parts[k].h, p->nodeID); */
Pedro Gonnet's avatar
Pedro Gonnet committed
885

886
      /* Re-link the gparts. */
887
888
      for (int kk = 0; kk < prox->nr_gparts_in; kk++) {
        struct gpart *gp = &s->gparts[offset_gparts + count_gparts + kk];
889
890

        if (gp->type == swift_type_gas) {
891
          struct part *p =
892
              &s->parts[offset_parts + count_parts - gp->id_or_neg_offset];