main.c 33.2 KB
Newer Older
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)
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
 *
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
 *
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
 *
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
 *
23
24
 ******************************************************************************/

Pedro Gonnet's avatar
Pedro Gonnet committed
25
26
/* Config parameters. */
#include "../config.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
27
28

/* Some standard headers. */
29
#include <errno.h>
30
#include <fenv.h>
31
#include <libgen.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
32
33
34
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
35
#include <sys/stat.h>
36
#include <unistd.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
37

38
39
/* MPI headers. */
#ifdef WITH_MPI
40
#include <mpi.h>
41
42
#endif

Pedro Gonnet's avatar
Pedro Gonnet committed
43
/* Local headers. */
44
#include "swift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
45

46
47
/* Engine policy flags. */
#ifndef ENGINE_POLICY
48
#define ENGINE_POLICY engine_policy_none
49
50
#endif

James Willis's avatar
James Willis committed
51
52
53
/* Global profiler. */
struct profiler prof;

54
55
56
/**
 * @brief Help messages for the command line parameters.
 */
57
58
void print_help_message() {

Matthieu Schaller's avatar
Matthieu Schaller committed
59
  printf("\nUsage: swift [OPTION]... PARAMFILE\n");
60
  printf("       swift_mpi [OPTION]... PARAMFILE\n\n");
61

62
  printf("Valid options are:\n");
63
64
65
66
  printf("  %2s %14s %s\n", "-a", "", "Pin runners using processor affinity.");
  printf("  %2s %14s %s\n", "-c", "",
         "Run with cosmological time integration.");
  printf("  %2s %14s %s\n", "-C", "", "Run with cooling.");
67
  printf(
68
      "  %2s %14s %s\n", "-d", "",
69
70
      "Dry run. Read the parameter file, allocate memory but does not read ");
  printf(
71
      "  %2s %14s %s\n", "", "",
72
      "the particles from ICs and exit before the start of time integration.");
73
  printf("  %2s %14s %s\n", "", "",
74
75
         "Allows user to check validy of parameter and IC files as well as "
         "memory limits.");
76
  printf("  %2s %14s %s\n", "-D", "",
77
78
         "Always drift all particles even the ones far from active particles. "
         "This emulates");
79
  printf("  %2s %14s %s\n", "", "",
80
         "Gadget-[23] and GIZMO's default behaviours.");
81
  printf("  %2s %14s %s\n", "-e", "",
82
         "Enable floating-point exceptions (debugging mode).");
83
  printf("  %2s %14s %s\n", "-f", "{int}",
84
         "Overwrite the CPU frequency (Hz) to be used for time measurements.");
85
  printf("  %2s %14s %s\n", "-g", "",
86
         "Run with an external gravitational potential.");
87
88
  printf("  %2s %14s %s\n", "-G", "", "Run with self-gravity.");
  printf("  %2s %14s %s\n", "-M", "",
89
         "Reconstruct the multipoles every time-step.");
90
  printf("  %2s %14s %s\n", "-n", "{int}",
91
92
         "Execute a fixed number of time steps. When unset use the time_end "
         "parameter to stop.");
93
94
95
  printf("  %2s %14s %s\n", "-P", "{sec:par:val}",
         "Set parameter value and overwrites values read from the parameters "
         "file. Can be used more than once.");
96
  printf("  %2s %14s %s\n", "-r", "", "Continue using restart files.");
97
98
99
  printf("  %2s %14s %s\n", "-s", "", "Run with hydrodynamics.");
  printf("  %2s %14s %s\n", "-S", "", "Run with stars.");
  printf("  %2s %14s %s\n", "-t", "{int}",
100
101
         "The number of threads to use on each MPI rank. Defaults to 1 if not "
         "specified.");
102
103
104
105
106
  printf("  %2s %14s %s\n", "-T", "", "Print timers every time-step.");
  printf("  %2s %14s %s\n", "-v", "[12]", "Increase the level of verbosity:");
  printf("  %2s %14s %s\n", "", "", "1: MPI-rank 0 writes,");
  printf("  %2s %14s %s\n", "", "", "2: All MPI-ranks write.");
  printf("  %2s %14s %s\n", "-y", "{int}",
107
         "Time-step frequency at which task graphs are dumped.");
108
109
  printf("  %2s %14s %s\n", "-Y", "{int}",
         "Time-step frequency at which threadpool tasks are dumped.");
110
  printf("  %2s %14s %s\n", "-h", "", "Print this help message and exit.");
111
112
113
  printf(
      "\nSee the file parameter_example.yml for an example of "
      "parameter file.\n");
114
115
}

Pedro Gonnet's avatar
Pedro Gonnet committed
116
117
118
119
/**
 * @brief Main routine that loads a few particles and generates some output.
 *
 */
120
121
int main(int argc, char *argv[]) {

122
  struct clocks_time tic, toc;
123
  struct engine e;
124

125
  int nr_nodes = 1, myrank = 0;
126

127
#ifdef WITH_MPI
128
  /* Start by initializing MPI. */
129
  int res = 0, prov = 0;
130
131
  if ((res = MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &prov)) !=
      MPI_SUCCESS)
132
133
    error("Call to MPI_Init failed with error %i.", res);
  if (prov != MPI_THREAD_MULTIPLE)
134
    error(
135
136
        "MPI does not provide the level of threading"
        " required (MPI_THREAD_MULTIPLE).");
137
  if ((res = MPI_Comm_size(MPI_COMM_WORLD, &nr_nodes)) != MPI_SUCCESS)
138
139
140
    error("MPI_Comm_size failed with error %i.", res);
  if ((res = MPI_Comm_rank(MPI_COMM_WORLD, &myrank)) != MPI_SUCCESS)
    error("Call to MPI_Comm_rank failed with error %i.", res);
141
142
143
144

  /* Make sure messages are stamped with the correct rank. */
  engine_rank = myrank;

145
146
147
  if ((res = MPI_Comm_set_errhandler(MPI_COMM_WORLD, MPI_ERRORS_RETURN)) !=
      MPI_SUCCESS)
    error("Call to MPI_Comm_set_errhandler failed with error %i.", res);
148
149
  if (myrank == 0)
    printf("[0000] [00000.0] main: MPI is up and running with %i node(s).\n\n",
Peter W. Draper's avatar
Peter W. Draper committed
150
           nr_nodes);
151
152
153
154
  if (nr_nodes == 1) {
    message("WARNING: you are running with one MPI rank.");
    message("WARNING: you should use the non-MPI version of this program.");
  }
155
  fflush(stdout);
156

157
#endif
158

159
160
  /* Welcome to SWIFT, you made the right choice */
  if (myrank == 0) greetings();
161

162
  int with_aff = 0;
163
  int dry_run = 0;
164
  int dump_tasks = 0;
165
  int dump_threadpool = 0;
166
  int nsteps = -2;
167
  int restart = 0;
168
169
  int with_cosmology = 0;
  int with_external_gravity = 0;
Tom Theuns's avatar
Tom Theuns committed
170
  int with_sourceterms = 0;
171
  int with_cooling = 0;
172
173
  int with_self_gravity = 0;
  int with_hydro = 0;
174
  int with_stars = 0;
175
  int with_fp_exceptions = 0;
176
  int with_drift_all = 0;
177
  int with_mpole_reconstruction = 0;
178
  int verbose = 0;
179
  int nr_threads = 1;
180
  int with_verbose_timers = 0;
181
182
  int nparams = 0;
  char *cmdparams[PARSER_MAX_NO_OF_PARAMS];
183
  char paramFileName[200] = "";
184
  char restartfile[200] = "";
185
186
187
  unsigned long long cpufreq = 0;

  /* Parse the parameters */
188
  int c;
189
  while ((c = getopt(argc, argv, "acCdDef:FgGhMn:P:rsSt:Tv:y:Y:")) != -1)
190
    switch (c) {
191
      case 'a':
192
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA)
193
        with_aff = 1;
194
195
196
#else
        error("Need NUMA support for thread affinity");
#endif
197
        break;
198
      case 'c':
199
        with_cosmology = 1;
200
        break;
201
202
203
      case 'C':
        with_cooling = 1;
        break;
204
      case 'd':
205
        dry_run = 1;
206
        break;
207
208
209
      case 'D':
        with_drift_all = 1;
        break;
210
      case 'e':
211
#ifdef HAVE_FE_ENABLE_EXCEPT
212
        with_fp_exceptions = 1;
213
214
215
#else
        error("Need support for floating point exception on this platform");
#endif
216
        break;
217
      case 'f':
218
219
220
        if (sscanf(optarg, "%llu", &cpufreq) != 1) {
          if (myrank == 0) printf("Error parsing CPU frequency (-f).\n");
          if (myrank == 0) print_help_message();
221
          return 1;
222
        }
223
        break;
224
225
226
      case 'F':
        with_sourceterms = 1;
        break;
227
      case 'g':
228
        with_external_gravity = 1;
229
        break;
230
231
      case 'G':
        with_self_gravity = 1;
Tom Theuns's avatar
Tom Theuns committed
232
        break;
233
      case 'h':
234
        if (myrank == 0) print_help_message();
235
        return 0;
236
237
238
      case 'M':
        with_mpole_reconstruction = 1;
        break;
239
240
241
242
243
244
245
      case 'n':
        if (sscanf(optarg, "%d", &nsteps) != 1) {
          if (myrank == 0) printf("Error parsing fixed number of steps.\n");
          if (myrank == 0) print_help_message();
          return 1;
        }
        break;
246
247
248
249
      case 'P':
        cmdparams[nparams] = optarg;
        nparams++;
        break;
250
251
252
      case 'r':
        restart = 1;
        break;
253
      case 's':
254
        with_hydro = 1;
255
        break;
256
257
258
      case 'S':
        with_stars = 1;
        break;
259
260
261
262
263
264
265
266
      case 't':
        if (sscanf(optarg, "%d", &nr_threads) != 1) {
          if (myrank == 0)
            printf("Error parsing the number of threads (-t).\n");
          if (myrank == 0) print_help_message();
          return 1;
        }
        break;
267
268
269
      case 'T':
        with_verbose_timers = 1;
        break;
270
      case 'v':
271
272
273
        if (sscanf(optarg, "%d", &verbose) != 1) {
          if (myrank == 0) printf("Error parsing verbosity level (-v).\n");
          if (myrank == 0) print_help_message();
274
          return 1;
275
        }
276
        break;
277
      case 'y':
278
279
280
        if (sscanf(optarg, "%d", &dump_tasks) != 1) {
          if (myrank == 0) printf("Error parsing dump_tasks (-y). \n");
          if (myrank == 0) print_help_message();
281
          return 1;
282
        }
283
284
285
286
287
288
#ifndef SWIFT_DEBUG_TASKS
        if (dump_tasks) {
          error(
              "Task dumping is only possible if SWIFT was configured with the "
              "--enable-task-debugging option.");
        }
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
#endif
        break;
      case 'Y':
        if (sscanf(optarg, "%d", &dump_threadpool) != 1) {
          if (myrank == 0) printf("Error parsing dump_threadpool (-Y). \n");
          if (myrank == 0) print_help_message();
          return 1;
        }
#ifndef SWIFT_DEBUG_THREADPOOL
        if (dump_threadpool) {
          error(
              "Threadpool dumping is only possible if SWIFT was configured "
              "with the "
              "--enable-threadpool-debugging option.");
        }
304
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
305
306
307
308
309
        break;
      case '?':
        if (myrank == 0) print_help_message();
        return 1;
        break;
310
    }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
311
312
313
314
315
  if (optind == argc - 1) {
    if (!strcpy(paramFileName, argv[optind++]))
      error("Error reading parameter file name.");
  } else if (optind > argc - 1) {
    if (myrank == 0) printf("Error: A parameter file name must be provided\n");
316
317
    if (myrank == 0) print_help_message();
    return 1;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
318
319
320
321
322
323
324
325
326
327
328
  } else {
    if (myrank == 0) printf("Error: Too many parameters given\n");
    if (myrank == 0) print_help_message();
    return 1;
  }
  if (!with_self_gravity && !with_hydro && !with_external_gravity) {
    if (myrank == 0)
      printf("Error: At least one of -s, -g or -G must be chosen.\n");
    if (myrank == 0) print_help_message();
    return 1;
  }
329
330
331
  if (with_stars && !with_external_gravity && !with_self_gravity) {
    if (myrank == 0)
      printf(
332
          "Error: Cannot process stars without gravity, -g or -G must be "
333
334
335
336
          "chosen.\n");
    if (myrank == 0) print_help_message();
    return 1;
  }
337

338
339
/* Let's pin the main thread, now we know if affinity will be used. */
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
340
341
  if (with_aff &&
      ((ENGINE_POLICY)&engine_policy_setaffinity) == engine_policy_setaffinity)
342
343
344
    engine_pin();
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
345
346
  /* Genesis 1.1: And then, there was time ! */
  clocks_set_cpufreq(cpufreq);
347

348
349
350
  /* How vocal are we ? */
  const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
351
352
353
  if (myrank == 0 && dry_run)
    message(
        "Executing a dry run. No i/o or time integration will be performed.");
354

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
355
356
357
358
359
  /* Report CPU frequency.*/
  cpufreq = clocks_get_cpufreq();
  if (myrank == 0) {
    message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
  }
360

Matthieu Schaller's avatar
Matthieu Schaller committed
361
/* Report host name(s). */
362
#ifdef WITH_MPI
363
  if (talking) {
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
364
365
    message("Rank %d running on: %s", myrank, hostname());
  }
366
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
367
  message("Running on: %s", hostname());
368
369
#endif

370
371
/* Do we have debugging checks ? */
#ifdef SWIFT_DEBUG_CHECKS
372
373
  if (myrank == 0)
    message("WARNING: Debugging checks activated. Code will be slower !");
374
375
#endif

376
377
378
379
380
381
382
/* Do we have debugging checks ? */
#ifdef SWIFT_USE_NAIVE_INTERACTIONS
  if (myrank == 0)
    message(
        "WARNING: Naive cell interactions activated. Code will be slower !");
#endif

383
384
385
386
387
388
389
390
391
/* Do we have gravity accuracy checks ? */
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
  if (myrank == 0)
    message(
        "WARNING: Checking 1/%d of all gpart for gravity accuracy. Code will "
        "be slower !",
        SWIFT_GRAVITY_FORCE_CHECKS);
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
392
393
  /* Do we choke on FP-exceptions ? */
  if (with_fp_exceptions) {
394
#ifdef HAVE_FE_ENABLE_EXCEPT
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
395
    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
396
#endif
397
398
    if (myrank == 0)
      message("WARNING: Floating point exceptions will be reported.");
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
399
  }
400

401
/* Do we have slow barriers? */
402
#ifndef HAVE_PTHREAD_BARRIERS
403
  if (myrank == 0)
404
    message("WARNING: Non-optimal thread barriers are being used.");
405
406
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
407
408
  /* How large are the parts? */
  if (myrank == 0) {
409
410
411
412
413
414
415
416
    message("sizeof(part)        is %4zi bytes.", sizeof(struct part));
    message("sizeof(xpart)       is %4zi bytes.", sizeof(struct xpart));
    message("sizeof(spart)       is %4zi bytes.", sizeof(struct spart));
    message("sizeof(gpart)       is %4zi bytes.", sizeof(struct gpart));
    message("sizeof(multipole)   is %4zi bytes.", sizeof(struct multipole));
    message("sizeof(grav_tensor) is %4zi bytes.", sizeof(struct grav_tensor));
    message("sizeof(task)        is %4zi bytes.", sizeof(struct task));
    message("sizeof(cell)        is %4zi bytes.", sizeof(struct cell));
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
417
  }
418

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
419
420
421
422
423
424
  /* Read the parameter file */
  struct swift_params *params = malloc(sizeof(struct swift_params));
  if (params == NULL) error("Error allocating memory for the parameter file.");
  if (myrank == 0) {
    message("Reading runtime parameters from file '%s'", paramFileName);
    parser_read_file(paramFileName, params);
425
426

    /* Handle any command-line overrides. */
427
428
429
430
    if (nparams > 0) {
      message(
          "Overwriting values read from the YAML file with command-line "
          "values.");
431
      for (int k = 0; k < nparams; k++) parser_set_param(params, cmdparams[k]);
432
    }
433
434

    /* And dump the parameters as used. */
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
435
436
437
    // parser_print_params(&params);
    parser_write_params_to_file(params, "used_parameters.yml");
  }
438
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
439
440
  /* Broadcast the parameter file */
  MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
441
#endif
442

443
  /* Check that we can write the snapshots by testing if the output
444
   * directory exists and is searchable and writable. */
445
446
447
  char basename[PARSER_MAX_LINE_SIZE];
  parser_get_param_string(params, "Snapshots:basename", basename);
  const char *dirp = dirname(basename);
448
449
  if (access(dirp, W_OK | X_OK) != 0) {
    error("Cannot write snapshots in directory %s (%s)", dirp, strerror(errno));
450
451
  }

452
  /* Prepare the domain decomposition scheme */
453
  struct repartition reparttype;
454
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
455
456
457
458
459
460
461
462
463
464
  struct partition initial_partition;
  partition_init(&initial_partition, &reparttype, params, nr_nodes);

  /* Let's report what we did */
  if (myrank == 0) {
    message("Using initial partition %s",
            initial_partition_name[initial_partition.type]);
    if (initial_partition.type == INITPART_GRID)
      message("grid set to [ %i %i %i ].", initial_partition.grid[0],
              initial_partition.grid[1], initial_partition.grid[2]);
465
    message("Using %s repartitioning", repartition_name[reparttype.type]);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
466
  }
467
#endif
468

469
470
471
472
473
  /* Common variables for restart and IC sections. */
  int clean_h_values = 0;
  int flag_entropy_ICs = 0;

  /* Work out where we will read and write restart files. */
474
  char restartdir[PARSER_MAX_LINE_SIZE];
475
  parser_get_opt_param_string(params, "Restarts:subdir", restartdir, "restart");
476
477
478

  /* The directory must exist. */
  if (access(restartdir, W_OK | X_OK) != 0) {
479
480
481
482
483
484
    if (restart) {
      error("Cannot restart as no restart subdirectory: %s (%s)", restartdir,
            strerror(errno));
    } else {
      if (mkdir(restartdir, 0777) != 0)
        error("Failed to create restart directory: %s (%s)", restartdir,
485
              strerror(errno));
486
    }
487
  }
488

489
490
491
492
  /* Basename for any restart files. */
  char restartname[PARSER_MAX_LINE_SIZE];
  parser_get_opt_param_string(params, "Restarts:basename", restartname,
                              "swift");
493

494
  if (restart) {
495

496
497
    /* Attempting a restart. */
    message("Restarting SWIFT");
498
499
500
501
502
    char **restartfiles = NULL;
    int nrestartfiles = 0;

    if (myrank == 0) {

503
      /* Locate the restart files. */
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
      restartfiles = restart_locate(restartdir, restartname, &nrestartfiles);
      if (nrestartfiles == 0)
        error("Failed to locate any restart files in %s", restartdir);

      /* We need one file per rank. */
      if (nrestartfiles != nr_nodes)
        error("Incorrect number of restart files, expected %d found %d",
              nr_nodes, nrestartfiles);

      if (verbose > 0)
        for (int i = 0; i < nrestartfiles; i++)
          message("found restart file: %s", restartfiles[i]);
    }

#ifdef WITH_MPI
519
    /* Distribute the restart files, need one for each rank. */
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
    if (myrank == 0) {

      for (int i = 1; i < nr_nodes; i++) {
        strcpy(restartfile, restartfiles[i]);
        MPI_Send(restartfile, 200, MPI_BYTE, i, 0, MPI_COMM_WORLD);
      }

      /* Keep local file. */
      strcpy(restartfile, restartfiles[0]);

      /* Finished with the list. */
      restart_locate_free(nrestartfiles, restartfiles);

    } else {
      MPI_Recv(restartfile, 200, MPI_BYTE, 0, 0, MPI_COMM_WORLD,
               MPI_STATUS_IGNORE);
    }
537
    if (verbose > 1) message("local restart file = %s", restartfile);
538
#else
539

540
541
542
543
    /* Just one restart file. */
    strcpy(restartfile, restartfiles[0]);
#endif

544
545
    /* Now read it. */
    restart_read(&e, restartfile);
546
547
548

    /* And initialize the engine with the space and policies. */
    if (myrank == 0) clocks_gettime(&tic);
549
    engine_config(1, &e, nr_nodes, myrank, nr_threads, with_aff, talking);
550
551
    if (myrank == 0) {
      clocks_gettime(&toc);
552
      message("engine_config took %.3f %s.", clocks_diff(&tic, &toc),
553
554
555
556
              clocks_getunit());
      fflush(stdout);
    }

557
558
559
560
561
    /* Check if we are already done when given steps on the command-line. */
    if (e.step >= nsteps && nsteps > 0)
      error("Not restarting, already completed %d steps (of out %d)", e.step,
            nsteps);

562
  } else {
563

564
565
566
567
568
569
570
571
572
573
574
575
576
577
    /* Reading ICs. */
    /* Initialize unit system and constants */
    struct unit_system us;
    struct phys_const prog_const;
    units_init(&us, params, "InternalUnitSystem");
    phys_const_init(&us, &prog_const);
    if (myrank == 0 && verbose > 0) {
      message("Internal unit system: U_M = %e g.", us.UnitMass_in_cgs);
      message("Internal unit system: U_L = %e cm.", us.UnitLength_in_cgs);
      message("Internal unit system: U_t = %e s.", us.UnitTime_in_cgs);
      message("Internal unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
      message("Internal unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
      phys_const_print(&prog_const);
    }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
578

579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
    /* Initialise the hydro properties */
    struct hydro_props hydro_properties;
    if (with_hydro) hydro_props_init(&hydro_properties, params);
    if (with_hydro) eos_init(&eos, params);

    /* Initialise the gravity properties */
    struct gravity_props gravity_properties;
    if (with_self_gravity) gravity_props_init(&gravity_properties, params);

    /* Read particles and space information from (GADGET) ICs */
    char ICfileName[200] = "";
    parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
    const int replicate =
        parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
    clean_h_values =
        parser_get_opt_param_int(params, "InitialConditions:cleanup_h", 0);
    if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
    fflush(stdout);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
597

598
599
600
601
602
603
604
605
    /* Get ready to read particles of all kinds */
    struct part *parts = NULL;
    struct gpart *gparts = NULL;
    struct spart *sparts = NULL;
    size_t Ngas = 0, Ngpart = 0, Nspart = 0;
    double dim[3] = {0., 0., 0.};
    int periodic = 0;
    if (myrank == 0) clocks_gettime(&tic);
606
607
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
608
609
610
611
612
613
614
    read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
                     &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                     (with_external_gravity || with_self_gravity), with_stars,
                     myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL,
                     nr_threads, dry_run);
#else
    read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
Matthieu Schaller's avatar
Matthieu Schaller committed
615
616
                   &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                   (with_external_gravity || with_self_gravity), with_stars,
617
618
                   myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
                   dry_run);
619
620
#endif
#else
621
622
623
624
    read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
                   &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                   (with_external_gravity || with_self_gravity), with_stars,
                   nr_threads, dry_run);
625
#endif
626
627
628
629
630
631
    if (myrank == 0) {
      clocks_gettime(&toc);
      message("Reading initial conditions took %.3f %s.",
              clocks_diff(&tic, &toc), clocks_getunit());
      fflush(stdout);
    }
632

633
#ifdef SWIFT_DEBUG_CHECKS
634
635
636
637
638
639
640
641
642
    /* Check once and for all that we don't have unwanted links */
    if (!with_stars) {
      for (size_t k = 0; k < Ngpart; ++k)
        if (gparts[k].type == swift_type_star) error("Linking problem");
    }
    if (!with_hydro) {
      for (size_t k = 0; k < Ngpart; ++k)
        if (gparts[k].type == swift_type_gas) error("Linking problem");
    }
643
#endif
644

645
646
    /* Get the total number of particles across all nodes. */
    long long N_total[3] = {0, 0, 0};
647
#if defined(WITH_MPI)
648
649
650
    long long N_long[3] = {Ngas, Ngpart, Nspart};
    MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
                  MPI_COMM_WORLD);
651
#else
652
653
654
    N_total[0] = Ngas;
    N_total[1] = Ngpart;
    N_total[2] = Nspart;
655
#endif
656
657
658
659
660
    if (myrank == 0)
      message(
          "Read %lld gas particles, %lld star particles and "
          "%lld gparts from the ICs.",
          N_total[0], N_total[2], N_total[1]);
661

662
663
664
665
666
    /* Initialize the space with these data. */
    if (myrank == 0) clocks_gettime(&tic);
    struct space s;
    space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
               periodic, replicate, with_self_gravity, talking, dry_run);
667

668
669
670
671
672
673
    if (myrank == 0) {
      clocks_gettime(&toc);
      message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
              clocks_getunit());
      fflush(stdout);
    }
674

675
676
677
678
679
680
681
682
683
684
685
686
687
    /* Say a few nice things about the space we just created. */
    if (myrank == 0) {
      message("space dimensions are [ %.3f %.3f %.3f ].", s.dim[0], s.dim[1],
              s.dim[2]);
      message("space %s periodic.", s.periodic ? "is" : "isn't");
      message("highest-level cell dimensions are [ %i %i %i ].", s.cdim[0],
              s.cdim[1], s.cdim[2]);
      message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
      message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
      message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
      message("maximum depth is %d.", s.maxdepth);
      fflush(stdout);
    }
688

689
690
691
692
693
694
    /* Verify that each particle is in it's proper cell. */
    if (talking && !dry_run) {
      int icount = 0;
      space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
      message("map_cellcheck picked up %i parts.", icount);
    }
695

696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
    /* Verify the maximal depth of cells. */
    if (talking && !dry_run) {
      int data[2] = {s.maxdepth, 0};
      space_map_cells_pre(&s, 0, &map_maxdepth, data);
      message("nr of cells at depth %i is %i.", data[0], data[1]);
    }

    /* Initialise the external potential properties */
    struct external_potential potential;
    if (with_external_gravity)
      potential_init(params, &prog_const, &us, &s, &potential);
    if (with_external_gravity && myrank == 0) potential_print(&potential);

    /* Initialise the cooling function properties */
    struct cooling_function_data cooling_func;
    if (with_cooling) cooling_init(params, &us, &prog_const, &cooling_func);
    if (with_cooling && myrank == 0) cooling_print(&cooling_func);

    /* Initialise the feedback properties */
    struct sourceterms sourceterms;
    if (with_sourceterms) sourceterms_init(params, &us, &sourceterms);
    if (with_sourceterms && myrank == 0) sourceterms_print(&sourceterms);

    /* Construct the engine policy */
    int engine_policies = ENGINE_POLICY | engine_policy_steal;
    if (with_drift_all) engine_policies |= engine_policy_drift_all;
    if (with_mpole_reconstruction)
      engine_policies |= engine_policy_reconstruct_mpoles;
    if (with_hydro) engine_policies |= engine_policy_hydro;
    if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
    if (with_external_gravity)
      engine_policies |= engine_policy_external_gravity;
    if (with_cosmology) engine_policies |= engine_policy_cosmology;
    if (with_cooling) engine_policies |= engine_policy_cooling;
    if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
    if (with_stars) engine_policies |= engine_policy_stars;

    /* Initialize the engine with the space and policies. */
    if (myrank == 0) clocks_gettime(&tic);
735
736
737
738
    engine_init(&e, &s, params, N_total[0], N_total[1], engine_policies,
                talking, &reparttype, &us, &prog_const, &hydro_properties,
                &gravity_properties, &potential, &cooling_func, &sourceterms);
    engine_config(0, &e, nr_nodes, myrank, nr_threads, with_aff, talking);
739
740
741
742
743
744
    if (myrank == 0) {
      clocks_gettime(&toc);
      message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
              clocks_getunit());
      fflush(stdout);
    }
745

746
#if defined(WITH_MPI)
747
748
749
750
751
    N_long[0] = s.nr_parts;
    N_long[1] = s.nr_gparts;
    N_long[2] = s.nr_sparts;
    MPI_Reduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM, 0,
               MPI_COMM_WORLD);
752
#else
753
754
755
    N_total[0] = s.nr_parts;
    N_total[1] = s.nr_gparts;
    N_total[2] = s.nr_sparts;
756
757
#endif

758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
    /* Get some info to the user. */
    if (myrank == 0) {
      long long N_DM = N_total[1] - N_total[2] - N_total[0];
      message(
          "Running on %lld gas particles, %lld star particles and %lld DM "
          "particles (%lld gravity particles)",
          N_total[0], N_total[2], N_total[1] > 0 ? N_DM : 0, N_total[1]);
      message(
          "from t=%.3e until t=%.3e with %d threads and %d queues "
          "(dt_min=%.3e, "
          "dt_max=%.3e)...",
          e.timeBegin, e.timeEnd, e.nr_threads, e.sched.nr_queues, e.dt_min,
          e.dt_max);
      fflush(stdout);
    }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
773
  }
774

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
775
776
  /* Time to say good-bye if this was not a serious run. */
  if (dry_run) {
777
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
778
779
    if ((res = MPI_Finalize()) != MPI_SUCCESS)
      error("call to MPI_Finalize failed with error %i.", res);
780
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
781
782
783
784
785
786
    if (myrank == 0)
      message("Time integration ready to start. End of dry-run.");
    engine_clean(&e);
    free(params);
    return 0;
  }
787

788
789
790
/* Initialise the table of Ewald corrections for the gravity checks */
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
  if (periodic) gravity_exact_force_ewald_init(e.s->dim[0]);
791
792
#endif

793
794
795
796
/* Init the runner history. */
#ifdef HIST
  for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
#endif
797

798
799
800
801
802
803
804
805
806
807
808
#ifdef WITH_MPI
  if (restart) {
    /* Set up the space and MPI */
    engine_makeproxies(&e);

  } else {
    /* Split the space. */
    engine_split(&e, &initial_partition);
    engine_redistribute(&e);
  }
#endif
809
  if (!restart) {
810
811
812
813
814
815
816
    /* Initialise the particles */
    engine_init_particles(&e, flag_entropy_ICs, clean_h_values);

    /* Write the state of the system before starting time integration. */
    engine_dump_snapshot(&e);
    engine_print_stats(&e);
  }
817

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
818
  /* Legend */
819
  if (myrank == 0)
820
    printf("# %6s %14s %14s %12s %12s %12s %16s [%s] %6s\n", "Step", "Time",
821
           "Time-step", "Updates", "g-Updates", "s-Updates", "Wall-clock time",
822
           clocks_getunit(), "Props");
823

824
  /* File for the timers */
825
  if (with_verbose_timers) timers_open_file(myrank);
826

827
828
829
830
  /* Create a name for restart file of this rank. */
  if (restart_genname(restartdir, restartname, e.nodeID, restartfile, 200) != 0)
    error("Failed to generate restart filename");

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
831
  /* Main simulation loop */
832
  for (int j = 0; !engine_is_done(&e) && e.step - 1 != nsteps; j++) {
833

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
834
    /* Reset timers */
835
    timers_reset_all();
836

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
837
838
    /* Take a step. */
    engine_step(&e);
839

840
    /* Print the timers. */
841
    if (with_verbose_timers) timers_print(e.step);
842

843
844
    /* Create restart files if required. Time/no. steps?*/
    restart_write(&e, restartfile);
845

846
#ifdef SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
847
848
    /* Dump the task data using the given frequency. */
    if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
849
850
#ifdef WITH_MPI

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
851
852
853
854
855
856
      /* Make sure output file is empty, only on one rank. */
      char dumpfile[30];
      snprintf(dumpfile, 30, "thread_info_MPI-step%d.dat", j + 1);
      FILE *file_thread;
      if (myrank == 0) {
        file_thread = fopen(dumpfile, "w");
857
        fclose(file_thread);
858
      }
859
      MPI_Barrier(MPI_COMM_WORLD);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
860
861
862
863
864
865
866
867
868
869
870
871

      for (int i = 0; i < nr_nodes; i++) {

        /* Rank 0 decides the index of writing node, this happens one-by-one. */
        int kk = i;
        MPI_Bcast(&kk, 1, MPI_INT, 0, MPI_COMM_WORLD);

        if (i == myrank) {

          /* Open file and position at end. */
          file_thread = fopen(dumpfile, "a");

Matthieu Schaller's avatar
Matthieu Schaller committed
872
873
874
          fprintf(file_thread, " %03i 0 0 0 0 %lli %lli %zi %zi %zi 0 0 %lli\n",
                  myrank, e.tic_step, e.toc_step, e.updates, e.g_updates,
                  e.s_updates, cpufreq);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
875
876
877
878
          int count = 0;
          for (int l = 0; l < e.sched.nr_tasks; l++) {
            if (!e.sched.tasks[l].implicit && e.sched.tasks[l].toc != 0) {
              fprintf(
Matthieu Schaller's avatar
Matthieu Schaller committed
879
880
881
                  file_thread,
                  " %03i %i %i %i %i %lli %lli %i %i %i %i %i %i\n", myrank,
                  e.sched.tasks[l].rid, e.sched.tasks[l].type,
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
882
883
884
885
886
887
888
889
890
891
                  e.sched.tasks[l].subtype, (e.sched.tasks[l].cj == NULL),
                  e.sched.tasks[l].tic, e.sched.tasks[l].toc,
                  (e.sched.tasks[l].ci != NULL) ? e.sched.tasks[l].ci->count
                                                : 0,
                  (e.sched.tasks[l].cj != NULL) ? e.sched.tasks[l].cj->count
                                                : 0,
                  (e.sched.tasks[l].ci != NULL) ? e.sched.tasks[l].ci->gcount
                                                : 0,
                  (e.sched.tasks[l].cj != NULL) ? e.sched.tasks[l].cj->gcount
                                                : 0,
892
                  e.sched.tasks[l].flags, e.sched.tasks[l].sid);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
893
894
895
896
897
898
899
900
901
902
            }
            fflush(stdout);
            count++;
          }
          fclose(file_thread);
        }

        /* And we wait for all to synchronize. */
        MPI_Barrier(MPI_COMM_WORLD);
      }
903

904
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
905
906
907
908
909
      char dumpfile[30];
      snprintf(dumpfile, 30, "thread_info-step%d.dat", j + 1);
      FILE *file_thread;
      file_thread = fopen(dumpfile, "w");
      /* Add some information to help with the plots */
910
911
912
      fprintf(file_thread, " %i %i %i %i %lli %lli %zi %zi %zi %i %lli\n", -2,
              -1, -1, 1, e.tic_step, e.toc_step, e.updates, e.g_updates,
              e.s_updates, 0, cpufreq);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
913
914
915
      for (int l = 0; l < e.sched.nr_tasks; l++) {
        if (!e.sched.tasks[l].implicit && e.sched.tasks[l].toc != 0) {
          fprintf(
916
              file_thread, " %i %i %i %i %lli %lli %i %i %i %i %i\n",
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
917
918
919
920
921
922
              e.sched.tasks[l].rid, e.sched.tasks[l].type,
              e.sched.tasks[l].subtype, (e.sched.tasks[l].cj == NULL),
              e.sched.tasks[l].tic, e.sched.tasks[l].toc,
              (e.sched.tasks[l].ci == NULL) ? 0 : e.sched.tasks[l].ci->count,
              (e.sched.tasks[l].cj == NULL) ? 0 : e.sched.tasks[l].cj->count,
              (e.sched.tasks[l].ci == NULL) ? 0 : e.sched.tasks[l].ci->gcount,
923
924
              (e.sched.tasks[l].cj == NULL) ? 0 : e.sched.tasks[l].cj->gcount,
              e.sched.tasks[l].sid);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
925
        }
926
      }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
927
      fclose(file_thread);
928
#endif  // WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
929
    }
930
#endif  // SWIFT_DEBUG_TASKS
931
932
933
934
935

#ifdef SWIFT_DEBUG_THREADPOOL
    /* Dump the task data using the given frequency. */
    if (dump_threadpool && (dump_threadpool == 1 || j % dump_threadpool == 1)) {
      char dumpfile[40];
936
#ifdef WITH_MPI
937
      snprintf(dumpfile, 40, "threadpool_info-rank%d-step%d.dat", engine_rank,
938
939
               j + 1);
#else
940
      snprintf(dumpfile, 40, "threadpool_info-step%d.dat", j + 1);
941
#endif  // WITH_MPI
942
943
944
945
946
      threadpool_dump_log(&e.threadpool, dumpfile, 1);
    } else {
      threadpool_reset_log(&e.threadpool);
    }
#endif  // SWIFT_DEBUG_THREADPOOL
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
947
  }
948
949

/* Print the values of the runner histogram. */
950
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
951
952
953
954
955
956
957
  printf("main: runner histogram data:\n");
  for (k = 0; k < runner_hist_N; k++)
    printf(" %e %e %e\n",
           runner_hist_a + k * (runner_hist_b - runner_hist_a) / runner_hist_N,
           runner_hist_a +
               (k + 1) * (runner_hist_b - runner_hist_a) / runner_hist_N,
           (double)runner_hist_bins[k]);
958
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
959

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
960
  /* Write final output. */
961
  engine_drift_all(&e);
962
  engine_print_stats(&e);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
963
  engine_dump_snapshot(&e);
964

965
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
966
967
  if ((res = MPI_Finalize()) != MPI_SUCCESS)
    error("call to MPI_Finalize failed with error %i.", res);
968
#endif
969

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
970
  /* Clean everything */
971
  if (with_verbose_timers) timers_close_file();
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
972
973
  engine_clean(&e);
  free(params);
974

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
975
976
  /* Say goodbye. */
  if (myrank == 0) message("done. Bye.");