engine.c 172 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
    if (e->policy & engine_policy_logger) {
460
      const uint32_t logger_flag = logger_pack_flags_and_data(
Loic Hausammann's avatar
Loic Hausammann committed
461
462
463
464
        logger_flag_mpi_exit, node_id);

      /* Log the particle when leaving a rank. */
      logger_log_part(e->logger, &s->parts[offset_parts + k],
465
                      &s->xparts[offset_parts + k],
Loic Hausammann's avatar
Loic Hausammann committed
466
467
468
469
                      logger_masks_all_part |
                      logger_mask_data[logger_special_flags].mask,
                      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
    if (e->policy & engine_policy_logger) {
509
      const uint32_t logger_flag = logger_pack_flags_and_data(
Loic Hausammann's avatar
Loic Hausammann committed
510
        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
      logger_log_spart(e->logger, &s->sparts[offset_sparts + k],
                       logger_masks_all_spart |
                       logger_mask_data[logger_special_flags].mask,
                       logger_flag);
    }
518
#endif
519
  }
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
  /* 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);
554
555

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

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

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

    /* Get the target node and proxy ID. */
569
570
571
572
    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];
573
    if (pid < 0) {
574
575
576
577
578
579
      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]);
580
581
582
583
584
585
    }

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

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

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

595
      const uint32_t logger_flag = logger_pack_flags_and_data(
Loic Hausammann's avatar
Loic Hausammann committed
596
         logger_flag_mpi_exit, node_id);
597
598

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

607
  /* Launch the proxies. */
608
609
  MPI_Request reqs_in[5 * engine_maxproxies];
  MPI_Request reqs_out[5 * engine_maxproxies];
610
  for (int k = 0; k < e->nr_proxies; k++) {
611
    proxy_parts_exchange_first(&e->proxies[k]);
612
613
614
615
616
    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. */
617
  for (int k = 0; k < e->nr_proxies; k++) {
618
    int pid = MPI_UNDEFINED;
619
620
    if (MPI_Waitany(e->nr_proxies, reqs_in, &pid, MPI_STATUS_IGNORE) !=
            MPI_SUCCESS ||
621
622
623
        pid == MPI_UNDEFINED)
      error("MPI_Waitany failed.");
    // message( "request from proxy %i has arrived." , pid );
624
    proxy_parts_exchange_second(&e->proxies[pid]);
625
  }
626

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

631
  /* Count the total number of incoming particles and make sure we have
632
     enough space to accommodate them. */
633
634
  int count_parts_in = 0;
  int count_gparts_in = 0;
635
  int count_sparts_in = 0;
636
  int count_bparts_in = 0;
637
638
639
  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;
640
    count_sparts_in += e->proxies[k].nr_sparts_in;
641
    count_bparts_in += e->proxies[k].nr_bparts_in;
642
643
  }
  if (e->verbose) {
644
645
646
647
648
    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);
649
  }
650
651

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

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

676
677
678
  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;
679
    if (swift_memalign("sparts", (void **)&sparts_new, spart_align,
680
681
682
                       sizeof(struct spart) * s->size_sparts) != 0)
      error("Failed to allocate new spart data.");
    memcpy(sparts_new, s->sparts, sizeof(struct spart) * offset_sparts);
683
    swift_free("sparts", s->sparts);
684
    s->sparts = sparts_new;
685
686

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

694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
  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;
      }
    }
  }
711

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

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

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

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

  /* Wait for each part array to come in and collect the new
     parts from the proxies. */
792
  int count_parts = 0, count_gparts = 0, count_sparts = 0, count_bparts = 0;
793
  for (int k = 0; k < nr_in; k++) {
794
    int err, pid;
795
    if ((err = MPI_Waitany(5 * e->nr_proxies, reqs_in, &pid,
796
                           MPI_STATUS_IGNORE)) != MPI_SUCCESS) {
797
798
799
800
801
802
      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;
803
804
    // message( "request from proxy %i has arrived." , pid / 5 );
    pid = 5 * (pid / 5);
Pedro Gonnet's avatar
Pedro Gonnet committed
805

806
    /* If all the requests for a given proxy have arrived... */
807
808
    if (reqs_in[pid + 0] == MPI_REQUEST_NULL &&
        reqs_in[pid + 1] == MPI_REQUEST_NULL &&
809
        reqs_in[pid + 2] == MPI_REQUEST_NULL &&
810
        reqs_in[pid + 3] == MPI_REQUEST_NULL &&
811
        reqs_in[pid + 4] == MPI_REQUEST_NULL) {
812
      /* Copy the particle data to the part/xpart/gpart arrays. */
813
      struct proxy *prox = &e->proxies[pid / 5];
814
815
816
817
818
819
      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);
820
821
      memcpy(&s->sparts[offset_sparts + count_sparts], prox->sparts_in,
             sizeof(struct spart) * prox->nr_sparts_in);
822
823
      memcpy(&s->bparts[offset_bparts + count_bparts], prox->bparts_in,
             sizeof(struct bpart) * prox->nr_bparts_in);
824
825

#ifdef WITH_LOGGER
Loic Hausammann's avatar
Loic Hausammann committed
826
      if (e->policy & engine_policy_logger) {
827
        const uint32_t flag = logger_pack_flags_and_data(logger_flag_mpi_enter,
Loic Hausammann's avatar
Loic Hausammann committed
828
829
830
831
832
833
834
835
                                              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 */
Loic Hausammann's avatar
Loic Hausammann committed
836
        logger_log_parts(e->logger, parts, xparts,
837
838
839
                         prox->nr_parts_in, logger_masks_all_part |
                         logger_mask_data[logger_special_flags].mask,
                         flag);
Loic Hausammann's avatar
Loic Hausammann committed
840
841

        /* Log the stellar particles */
842
843
844
845
        logger_log_sparts(e->logger, sparts, prox->nr_sparts_in,
                          logger_masks_all_spart |
                          logger_mask_data[logger_special_flags].mask,
                          flag);
Loic Hausammann's avatar
Loic Hausammann committed
846
847

        /* Log the gparts */
848
849
850
851
        logger_log_gparts(e->logger, gparts, prox->nr_gparts_in,
                          logger_masks_all_gpart |
                          logger_mask_data[logger_special_flags].mask,
                          flag);
Loic Hausammann's avatar
Loic Hausammann committed
852
853
854
855
856
857

        /* Log the bparts */
        if (prox->nr_bparts_in > 0) {
          error("TODO");
        }
      }
858
#endif
859