main.c 34.5 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 restart_file[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
475
  char restart_dir[PARSER_MAX_LINE_SIZE];
  parser_get_opt_param_string(params, "Restarts:subdir", restart_dir, "restart");
476
477

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

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

496
497
498
499
500
  /* How often to check for the stop file and dump restarts and exit the
   * application. */
  int restart_stop_steps = parser_get_opt_param_int(params, "Restarts:stop_steps", 100);

  /* If restarting, look for the restart files. */
501
  if (restart) {
502

503
    /* Attempting a restart. */
504
505
    char **restart_files = NULL;
    int restart_nfiles = 0;
506
507

    if (myrank == 0) {
508
      message("Restarting SWIFT");
509

510
      /* Locate the restart files. */
511
512
513
      restart_files = restart_locate(restart_dir, restart_name, &restart_nfiles);
      if (restart_nfiles == 0)
        error("Failed to locate any restart files in %s", restart_dir);
514
515

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

      if (verbose > 0)
521
522
        for (int i = 0; i < restart_nfiles; i++)
          message("found restart file: %s", restart_files[i]);
523
524
525
    }

#ifdef WITH_MPI
526
    /* Distribute the restart files, need one for each rank. */
527
528
529
    if (myrank == 0) {

      for (int i = 1; i < nr_nodes; i++) {
530
531
        strcpy(restart_file, restart_files[i]);
        MPI_Send(restart_file, 200, MPI_BYTE, i, 0, MPI_COMM_WORLD);
532
533
534
      }

      /* Keep local file. */
535
      strcpy(restart_file, restart_files[0]);
536
537

      /* Finished with the list. */
538
      restart_locate_free(restart_nfiles, restart_files);
539
540

    } else {
541
      MPI_Recv(restart_file, 200, MPI_BYTE, 0, 0, MPI_COMM_WORLD,
542
543
               MPI_STATUS_IGNORE);
    }
544
    if (verbose > 1) message("local restart file = %s", restart_file);
545
#else
546

547
    /* Just one restart file. */
548
    strcpy(restart_file, restart_files[0]);
549
550
#endif

551
    /* Now read it. */
552
    restart_read(&e, restart_file);
553
554
555

    /* And initialize the engine with the space and policies. */
    if (myrank == 0) clocks_gettime(&tic);
556
    engine_config(1, &e, params, nr_nodes, myrank, nr_threads, with_aff, talking,
557
                  restart_file);
558
559
    if (myrank == 0) {
      clocks_gettime(&toc);
560
      message("engine_config took %.3f %s.", clocks_diff(&tic, &toc),
561
562
563
564
              clocks_getunit());
      fflush(stdout);
    }

565
566
    /* Check if we are already done when given steps on the command-line. */
    if (e.step >= nsteps && nsteps > 0)
567
      error("Not restarting, already completed %d steps", e.step);
568

569
  } else {
570

571
    /* Not restarting so look for the ICs. */
572
573
574
575
576
577
578
579
580
581
582
583
584
    /* 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
585

586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
    /* 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
604

605
606
607
608
609
610
611
612
    /* 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);
613
614
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
615
616
617
618
619
620
621
    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
622
623
                   &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                   (with_external_gravity || with_self_gravity), with_stars,
624
625
                   myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
                   dry_run);
626
627
#endif
#else
628
629
630
631
    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);
632
#endif
633
634
635
636
637
638
    if (myrank == 0) {
      clocks_gettime(&toc);
      message("Reading initial conditions took %.3f %s.",
              clocks_diff(&tic, &toc), clocks_getunit());
      fflush(stdout);
    }
639

640
#ifdef SWIFT_DEBUG_CHECKS
641
642
643
644
645
646
647
648
649
    /* 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");
    }
650
#endif
651

652
653
    /* Get the total number of particles across all nodes. */
    long long N_total[3] = {0, 0, 0};
654
#if defined(WITH_MPI)
655
656
657
    long long N_long[3] = {Ngas, Ngpart, Nspart};
    MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
                  MPI_COMM_WORLD);
658
#else
659
660
661
    N_total[0] = Ngas;
    N_total[1] = Ngpart;
    N_total[2] = Nspart;
662
#endif
663

664
665
    if (myrank == 0)
      message(
666
667
          "Read %lld gas particles, %lld star particles and %lld gparts from the "
          "ICs.",
668
          N_total[0], N_total[2], N_total[1]);
669

670
671
672
673
674
    /* 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);
675

676
677
678
679
680
681
    if (myrank == 0) {
      clocks_gettime(&toc);
      message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
              clocks_getunit());
      fflush(stdout);
    }
682

683
    /* Also update the total counts (in case of changes due to replication) */
684
#if defined(WITH_MPI)
685
686
687
688
689
    N_long[0] = s.nr_parts;
    N_long[1] = s.nr_gparts;
    N_long[2] = s.nr_sparts;
    MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
                  MPI_COMM_WORLD);
690
#else
691
692
693
    N_total[0] = s.nr_parts;
    N_total[1] = s.nr_gparts;
    N_total[2] = s.nr_sparts;
694
695
#endif

696
697
698
699
700
701
702
703
704
705
706
707
708
    /* 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);
    }
709

710
711
712
713
714
715
    /* 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);
    }
716

717
718
719
720
721
722
723
724
725
726
727
    /* 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);
728
    if (myrank == 0) potential_print(&potential);
729
730
731
732

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

735
736
737
738
739
    /* Initialise the chemistry */
    struct chemistry_data chemistry;
    chemistry_init(params, &us, &prog_const, &chemistry);
    if (myrank == 0) chemistry_print(&chemistry);

740
741
742
    /* Initialise the feedback properties */
    struct sourceterms sourceterms;
    if (with_sourceterms) sourceterms_init(params, &us, &sourceterms);
Peter W. Draper's avatar
Peter W. Draper committed
743
    if (with_sourceterms && myrank == 0) sourceterms_print(&sourceterms);
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760

    /* 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);
761
762
    engine_init(&e, &s, params, N_total[0], N_total[1], engine_policies,
                talking, &reparttype, &us, &prog_const, &hydro_properties,
Peter W. Draper's avatar
Peter W. Draper committed
763
764
                &gravity_properties, &potential, &cooling_func, &chemistry,
                &sourceterms);
765
    engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff, talking,
766
                  restart_file);
767
768
769
770
771
772
    if (myrank == 0) {
      clocks_gettime(&toc);
      message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
              clocks_getunit());
      fflush(stdout);
    }
773

774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
    /* 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
789
  }
790

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
791
792
  /* Time to say good-bye if this was not a serious run. */
  if (dry_run) {
793
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
794
795
    if ((res = MPI_Finalize()) != MPI_SUCCESS)
      error("call to MPI_Finalize failed with error %i.", res);
796
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
797
798
799
800
801
802
    if (myrank == 0)
      message("Time integration ready to start. End of dry-run.");
    engine_clean(&e);
    free(params);
    return 0;
  }
803

804
805
806
/* 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]);
807
808
#endif

809
810
811
812
/* Init the runner history. */
#ifdef HIST
  for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
#endif
813

Peter W. Draper's avatar
Peter W. Draper committed
814
815
  if (!restart) {

816
817
818
819
820
#ifdef WITH_MPI
    /* Split the space. */
    engine_split(&e, &initial_partition);
    engine_redistribute(&e);
#endif
Peter W. Draper's avatar
Peter W. Draper committed
821

822
823
824
825
826
827
828
    /* 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);
  }
829

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
830
  /* Legend */
831
  if (myrank == 0)
832
    printf("# %6s %14s %14s %12s %12s %12s %16s [%s] %6s\n", "Step", "Time",
833
           "Time-step", "Updates", "g-Updates", "s-Updates", "Wall-clock time",
834
           clocks_getunit(), "Props");
835

836
  /* File for the timers */
837
  if (with_verbose_timers) timers_open_file(myrank);
838

839
  /* Create a name for restart file of this rank. */
840
  if (restart_genname(restart_dir, restart_name, e.nodeID, restart_file, 200) != 0)
841
842
    error("Failed to generate restart filename");

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
843
  /* Main simulation loop */
844
845
846
  /* ==================== */
  int force_stop = 0;
  for (int j = 0; !engine_is_done(&e) && e.step - 1 != nsteps && !force_stop; j++) {
847

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
848
    /* Reset timers */
849
    timers_reset_all();
850

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
851
852
    /* Take a step. */
    engine_step(&e);
853

854
    /* Print the timers. */
855
    if (with_verbose_timers) timers_print(e.step);
856

857
858
859
860
    /* Every so often allow the user to stop the application and dump the restart
     * files. */
    if (j % restart_stop_steps == 0) {
        force_stop = restart_stop_now(restart_dir, 0);
861
        if (myrank == 0 && force_stop)
862
863
864
865
            message("Forcing application exit, dumping restart files...");
    }

    /* Also if using nsteps to exit, will not have saved any restarts on exit, make
866
     * sure we do that (useful in testing only). */
867
    if (force_stop || (e.restart_onexit && e.step - 1 == nsteps)) engine_dump_restarts(&e, 0, 1);
868

869
#ifdef SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
870
871
    /* Dump the task data using the given frequency. */
    if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
872
873
#ifdef WITH_MPI

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
874
875
876
877
878
879
      /* 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");
880
        fclose(file_thread);
881
      }
882
      MPI_Barrier(MPI_COMM_WORLD);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
883
884
885
886
887
888
889
890
891
892
893
894

      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
895
896
897
          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
898
899
900
901
          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
902
903
904
                  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
905
906
907
908
909
910
911
912
913
914
                  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,
915
                  e.sched.tasks[l].flags, e.sched.tasks[l].sid);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
916
917
918
919
920
921
922
923
924
925
            }
            fflush(stdout);
            count++;
          }
          fclose(file_thread);
        }

        /* And we wait for all to synchronize. */
        MPI_Barrier(MPI_COMM_WORLD);
      }
Pedro Gonnet's avatar