main.c 42.4 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 "argparse.h"
45
#include "swift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
46

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

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

55
56
//  Usage string.
static const char *const swift_usage[] = {
57
58
59
60
61
62
63
  "swift [options] [[--] param-file]",
  "swift [options] param-file",
  "swift_mpi [options] [[--] param-file]",
  "swift_mpi [options] param-file",
  NULL,
};

64
65
66
67
68
69
70
71
72
73
74
75
76
// Function to handle multiple -P arguments.
struct cmdparams {
  const char *param[PARSER_MAX_NO_OF_PARAMS];
  int nparam;
};

static int handle_cmdparam(struct argparse *self,
                           const struct argparse_option *opt) {
  struct cmdparams *cmdps = (struct cmdparams *)opt->data;
  cmdps->param[cmdps->nparam] = *(char **)opt->value;
  cmdps->nparam++;
  return 1;
}
77

Pedro Gonnet's avatar
Pedro Gonnet committed
78
79
80
81
/**
 * @brief Main routine that loads a few particles and generates some output.
 *
 */
82
83
int main(int argc, char *argv[]) {

84
  struct clocks_time tic, toc;
85
  struct engine e;
86

87
88
  /* Structs used by the engine. Declare now to make sure these are always in
   * scope.  */
lhausamm's avatar
lhausamm committed
89
  struct chemistry_global_data chemistry;
90
  struct cooling_function_data cooling_func;
91
  struct cosmology cosmo;
92
  struct external_potential potential;
93
  struct pm_mesh mesh;
94
95
96
  struct gpart *gparts = NULL;
  struct gravity_props gravity_properties;
  struct hydro_props hydro_properties;
Loic Hausammann's avatar
Loic Hausammann committed
97
  struct stars_props stars_properties;
98
99
100
101
102
103
  struct part *parts = NULL;
  struct phys_const prog_const;
  struct sourceterms sourceterms;
  struct space s;
  struct spart *sparts = NULL;
  struct unit_system us;
104

105
  int nr_nodes = 1, myrank = 0;
106

107
#ifdef WITH_MPI
108
  /* Start by initializing MPI. */
109
  int res = 0, prov = 0;
110
111
  if ((res = MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &prov)) !=
      MPI_SUCCESS)
112
113
    error("Call to MPI_Init failed with error %i.", res);
  if (prov != MPI_THREAD_MULTIPLE)
114
    error(
115
116
        "MPI does not provide the level of threading"
        " required (MPI_THREAD_MULTIPLE).");
117
  if ((res = MPI_Comm_size(MPI_COMM_WORLD, &nr_nodes)) != MPI_SUCCESS)
118
119
120
    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);
121
122
123
124

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

125
126
127
  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);
128
129
  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
130
           nr_nodes);
131
132
133
134
  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.");
  }
135
  fflush(stdout);
136

137
#endif
138

139
140
  /* Welcome to SWIFT, you made the right choice */
  if (myrank == 0) greetings();
141

142
  int with_aff = 0;
143
  int dry_run = 0;
144
  int dump_tasks = 0;
145
  int dump_threadpool = 0;
146
  int nsteps = -2;
147
  int restart = 0;
148
149
  int with_cosmology = 0;
  int with_external_gravity = 0;
Tom Theuns's avatar
Tom Theuns committed
150
  int with_sourceterms = 0;
151
  int with_cooling = 0;
152
153
  int with_self_gravity = 0;
  int with_hydro = 0;
154
  int with_stars = 0;
Loic Hausammann's avatar
Loic Hausammann committed
155
  int with_feedback = 0;
156
  int with_fp_exceptions = 0;
157
  int with_drift_all = 0;
158
  int with_mpole_reconstruction = 0;
159
  int with_structure_finding = 0;
160
  int verbose = 0;
161
  int nr_threads = 1;
162
  int with_verbose_timers = 0;
163
  char *output_parameters_filename = NULL;
164
  char *cpufreqarg = NULL;
165
  char *param_filename = NULL;
166
  char restart_file[200] = "";
167
  unsigned long long cpufreq = 0;
168
169
170
171
  struct cmdparams cmdps;
  cmdps.nparam = 0;
  cmdps.param[0] = NULL;
  char *buffer = NULL;
172

173
174
175
  /* Parse the command-line parameters. */
  struct argparse_option options[] = {
      OPT_HELP(),
176

177
      OPT_GROUP("Simulation options"),
178
179
180
181
182
      OPT_BOOLEAN('b', "feedback", &with_feedback, "Run with stars feedback",
                  NULL, 0, 0),
      OPT_BOOLEAN('c', "cosmology", &with_cosmology,
                  "Run with cosmological time integration.", NULL, 0, 0),
      OPT_BOOLEAN('C', "cooling", &with_cooling, "Run with cooling", NULL, 0, 0),
183
      OPT_BOOLEAN('F', "sourceterms", &with_sourceterms, "", NULL, 0, 0),
184
185
186
187
188
189
190
191
192
193
194
195
      OPT_BOOLEAN('g', "external-gravity", &with_external_gravity,
                  "Run with an external gravitational potential.", NULL, 0, 0),
      OPT_BOOLEAN('G', "self-gravity", &with_self_gravity,
                  "Run with self-gravity.", NULL, 0, 0),
      OPT_BOOLEAN('M', "multipole-reconstruction", &with_mpole_reconstruction,
                  "Reconstruct the multipoles every time-step.", NULL, 0, 0),
      OPT_BOOLEAN('s', "hydrodynamics", &with_hydro, "Run with hydrodynamics.",
                  NULL, 0, 0),
      OPT_BOOLEAN('S', "stars", &with_stars, "Run with stars", NULL, 0, 0),
      OPT_BOOLEAN('x', "velociraptor", &with_structure_finding, 
                  "Run with structure finding", NULL, 0, 0),

196
      OPT_GROUP("Control options"),
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
      OPT_BOOLEAN('a', "affinity", &with_aff,
                  "Pin runners using processor affinity.", NULL, 0, 0),
      OPT_BOOLEAN('d', "dry-run", &dry_run, 
                  "Dry run. Read the parameter file, allocates memory but "
                  "does not read the particles from ICs.\t\nExits before the "
                  "start of time integration. Checks the validity of "
                  "parameters and IC files as well as memory limits.",
                  NULL, 0, 0),
      OPT_BOOLEAN('D', "drift-all", &with_drift_all, "Always drift all "
                  "particles even the ones far from active particles. "
                  "This emulates Gadget-[23] and GIZMO's default behaviours.",
                  NULL, 0, 0),
      OPT_BOOLEAN('e', "fpe", &with_fp_exceptions, "Enable floating-point"
                  "exceptions (debugging mode).", NULL, 0, 0),
      OPT_STRING('f', "cpu-frequency", &cpufreqarg, "Overwrite the CPU "
                 "frequency (Hz) to be used for time measurements.", NULL, 0, 0),
      OPT_INTEGER('n', "steps", &nsteps, "Execute a fixed number of time "
                  "steps. When unset use the time_end parameter to stop.",
                  NULL, 0, 0),
      OPT_STRING('o', "output-params", &output_parameters_filename, 
                 "Generate a default output parameter file.", NULL, 0, 0),
      OPT_STRING('P', "param", &buffer, "Set parameter value, "
                 "overwriting a value read from the parameter "
                 "file. Can be used more than once {sec:par:value}.",
                 handle_cmdparam, (intptr_t)&cmdps, 0),
      OPT_BOOLEAN('r', "restart", &restart, "Continue using restart files.",
                  NULL, 0, 0),
      OPT_INTEGER('t', "threads", &nr_threads, "The number of threads to "
                  "use on each MPI rank. Defaults to 1 if not specified.",
                  NULL, 0, 0),
      OPT_INTEGER('T', "timers", &with_verbose_timers,
               "Print timers every time-step.", NULL, 0, 0),
      OPT_INTEGER('v', "verbose", &verbose, "Run in verbose mode, in MPI "
                  "mode 2 outputs from all ranks.", NULL, 0, 0),
      OPT_INTEGER('y', "task-dumps", &dump_tasks, 
                  "Time-step frequency at which task graphs are dumped.",
                  NULL, 0, 0),
      OPT_INTEGER('Y', "threadpool-dumps", &dump_threadpool, 
                  "Time-step frequency at which threadpool tasks are dumped.",
                  NULL, 0, 0),
237
238
  };
  struct argparse argparse;
239
240
  argparse_init(&argparse, options, swift_usage, 0);
  argparse_describe(&argparse, "\nParameters:", NULL /* no epilog */);
241
242
243
244
  int nargs = argparse_parse(&argparse, argc, (const char **)argv);

  /* Need a parameter file. */
  if (nargs != 1) {
245
246
247
      if (myrank == 0) argparse_usage(&argparse);
      printf("\nError: no parameter file was supplied.\n");
      return 1;
248
249
250
251
252
  }
  param_filename = argv[0];

  /* Checks of options. */
#if ! defined(HAVE_SETAFFINITY) || !defined(HAVE_LIBNUMA)
253
254
255
256
  if (with_aff) {
      printf("Error: no NUMA support for thread affinity\n");
      return 1;
  }
257
#endif
258
259

#ifndef HAVE_FE_ENABLE_EXCEPT
260
261
262
263
  if (with_fp_exceptions) {
      printf("Error: no support for floating point exceptions\n");
      return 1;
  }
264
#endif
265
266

#ifndef HAVE_VELOCIRAPTOR
267
268
269
270
  if (with_structure_finding) {
      printf("Error: VELOCIraptor is not available\n");
      return 1;
  }
271
#endif
272

273
#ifndef SWIFT_DEBUG_TASKS
274
275
276
277
278
  if (dump_tasks) {
      printf("Error: task dumping is only possible if SWIFT was configured"
             " with the --enable-task-debugging option.\n");
      return 1;
  }
279
#endif
280

281
#ifndef SWIFT_DEBUG_THREADPOOL
282
283
284
285
286
  if (dump_threadpool) {
      printf("Error: threadpool dumping is only possible if SWIFT was "
             "configured with the --enable-threadpool-debugging option.\n");
      return 1;
  }
287
#endif
288

289
290
291
292
293
294
295
296
297
  /* The CPU frequency is a long long, so we need to parse that ourselves. */
  if (cpufreqarg != NULL) {
    if (sscanf(cpufreqarg, "%llu", &cpufreq) != 1) {
      if (myrank == 0) printf("Error parsing CPU frequency (%s).\n",
                              cpufreqarg);
      return 1;
    }
  }

298
  /* Write output parameter file */
299
  if (myrank == 0 && output_parameters_filename != NULL) {
300
    io_write_output_field_parameter(output_parameters_filename);
301
    printf("End of run.\n");
302
303
304
    return 0;
  }

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
305
  if (!with_self_gravity && !with_hydro && !with_external_gravity) {
306
307
308
309
310
      if (myrank == 0) {
          argparse_usage(&argparse);
          printf("\nError: At least one of -s, -g or -G must be chosen.\n");
      }
      return 1;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
311
  }
312
  if (with_stars && !with_external_gravity && !with_self_gravity) {
313
314
315
316
317
318
      if (myrank == 0) {
          argparse_usage(&argparse);
          printf("\nError: Cannot process stars without gravity, -g or -G "
                 "must be chosen.\n");
      }
      return 1;
319
  }
320

Loic Hausammann's avatar
Loic Hausammann committed
321
  if (!with_stars && with_feedback) {
322
323
324
325
326
327
      if (myrank == 0) {
          argparse_usage(&argparse);
          printf("\nError: Cannot process feedback without stars, -S must be "
                 "chosen.\n");
      }
      return 1;
Loic Hausammann's avatar
Loic Hausammann committed
328
329
  }

330
331
/* Let's pin the main thread, now we know if affinity will be used. */
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
332
333
  if (with_aff &&
      ((ENGINE_POLICY)&engine_policy_setaffinity) == engine_policy_setaffinity)
334
335
336
    engine_pin();
#endif

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

340
341
342
  /* How vocal are we ? */
  const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
347
348
349
350
351
  /* Report CPU frequency.*/
  cpufreq = clocks_get_cpufreq();
  if (myrank == 0) {
    message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
  }
352

Matthieu Schaller's avatar
Matthieu Schaller committed
353
/* Report host name(s). */
354
#ifdef WITH_MPI
355
  if (talking) {
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
356
357
    message("Rank %d running on: %s", myrank, hostname());
  }
358
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
359
  message("Running on: %s", hostname());
360
361
#endif

362
363
/* Do we have debugging checks ? */
#ifdef SWIFT_DEBUG_CHECKS
364
365
  if (myrank == 0)
    message("WARNING: Debugging checks activated. Code will be slower !");
366
367
#endif

368
369
370
371
372
373
374
/* Do we have debugging checks ? */
#ifdef SWIFT_USE_NAIVE_INTERACTIONS
  if (myrank == 0)
    message(
        "WARNING: Naive cell interactions activated. Code will be slower !");
#endif

375
376
377
378
379
380
381
382
383
/* 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
384
385
  /* Do we choke on FP-exceptions ? */
  if (with_fp_exceptions) {
386
#ifdef HAVE_FE_ENABLE_EXCEPT
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
387
    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
388
#endif
389
390
    if (myrank == 0)
      message("WARNING: Floating point exceptions will be reported.");
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
391
  }
392

393
/* Do we have slow barriers? */
394
#ifndef HAVE_PTHREAD_BARRIERS
395
  if (myrank == 0)
396
    message("WARNING: Non-optimal thread barriers are being used.");
397
398
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
399
400
  /* How large are the parts? */
  if (myrank == 0) {
401
402
403
404
405
406
407
408
    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
409
  }
410

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
411
  /* Read the parameter file */
Matthieu Schaller's avatar
Matthieu Schaller committed
412
413
  struct swift_params *params =
      (struct swift_params *)malloc(sizeof(struct swift_params));
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
414
415
  if (params == NULL) error("Error allocating memory for the parameter file.");
  if (myrank == 0) {
416
417
    message("Reading runtime parameters from file '%s'", param_filename);
    parser_read_file(param_filename, params);
418
419

    /* Handle any command-line overrides. */
420
    if (cmdps.nparam > 0) {
421
422
423
      message(
          "Overwriting values read from the YAML file with command-line "
          "values.");
424
      for (int k = 0; k < cmdps.nparam; k++) parser_set_param(params, cmdps.param[k]);
425
    }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
426
  }
427
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
428
429
  /* Broadcast the parameter file */
  MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
430
#endif
431

432
433
  /* Temporary early aborts for modes not supported over MPI. */
#ifdef WITH_MPI
Matthieu Schaller's avatar
Matthieu Schaller committed
434
435
  if (with_mpole_reconstruction && nr_nodes > 1)
    error("Cannot reconstruct m-poles every step over MPI (yet).");
436
437
#endif

438
439
440
441
#ifdef WITH_MPI
  if (with_feedback) error("Can't run with feedback over MPI (yet).");
#endif

442
#if defined(WITH_MPI) && defined(HAVE_VELOCIRAPTOR)
Matthieu Schaller's avatar
Matthieu Schaller committed
443
444
  if (with_structure_finding && nr_nodes > 1)
    error("VEOCIraptor not yet enabled over MPI.");
445
446
#endif

447
  /* Check that we can write the snapshots by testing if the output
448
   * directory exists and is searchable and writable. */
449
450
451
  char basename[PARSER_MAX_LINE_SIZE];
  parser_get_param_string(params, "Snapshots:basename", basename);
  const char *dirp = dirname(basename);
452
453
  if (access(dirp, W_OK | X_OK) != 0) {
    error("Cannot write snapshots in directory %s (%s)", dirp, strerror(errno));
454
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
455
456
457

  /* Check that we can write the structure finding catalogues by testing if the
   * output
458
   * directory exists and is searchable and writable. */
Matthieu Schaller's avatar
Matthieu Schaller committed
459
  if (with_structure_finding) {
460
    char stfbasename[PARSER_MAX_LINE_SIZE];
461
    parser_get_param_string(params, "StructureFinding:basename", stfbasename);
462
463
    const char *stfdirp = dirname(stfbasename);
    if (access(stfdirp, W_OK | X_OK) != 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
464
465
      error("Cannot write stf catalogues in directory %s (%s)", stfdirp,
            strerror(errno));
466
467
    }
  }
468

469
  /* Prepare the domain decomposition scheme */
470
  struct repartition reparttype;
471
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
472
473
474
475
476
  struct partition initial_partition;
  partition_init(&initial_partition, &reparttype, params, nr_nodes);

  /* Let's report what we did */
  if (myrank == 0) {
477
478
#if defined(HAVE_PARMETIS)
    if (reparttype.usemetis)
Matthieu Schaller's avatar
Matthieu Schaller committed
479
      message("Using METIS serial partitioning:");
480
    else
Matthieu Schaller's avatar
Matthieu Schaller committed
481
      message("Using ParMETIS partitioning:");
482
483
484
485
#else
    message("Using METIS serial partitioning:");
#endif
    message("  initial partitioning: %s",
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
486
487
            initial_partition_name[initial_partition.type]);
    if (initial_partition.type == INITPART_GRID)
488
      message("    grid set to [ %i %i %i ].", initial_partition.grid[0],
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
489
              initial_partition.grid[1], initial_partition.grid[2]);
490
    message("  repartitioning: %s", repartition_name[reparttype.type]);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
491
  }
492
#endif
493

494
  /* Common variables for restart and IC sections. */
495
  int clean_smoothing_length_values = 0;
496
497
498
  int flag_entropy_ICs = 0;

  /* Work out where we will read and write restart files. */
499
  char restart_dir[PARSER_MAX_LINE_SIZE];
Peter W. Draper's avatar
Peter W. Draper committed
500
501
  parser_get_opt_param_string(params, "Restarts:subdir", restart_dir,
                              "restart");
502
503

  /* The directory must exist. */
504
  if (myrank == 0) {
505
    if (access(restart_dir, W_OK | X_OK) != 0) {
506
      if (restart) {
507
        error("Cannot restart as no restart subdirectory: %s (%s)", restart_dir,
508
              strerror(errno));
509
      } else {
Peter W. Draper's avatar
Peter W. Draper committed
510
511
512
        if (mkdir(restart_dir, 0777) != 0)
          error("Failed to create restart directory: %s (%s)", restart_dir,
                strerror(errno));
513
      }
514
    }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
515
516
  }

517
  /* Basename for any restart files. */
518
519
  char restart_name[PARSER_MAX_LINE_SIZE];
  parser_get_opt_param_string(params, "Restarts:basename", restart_name,
520
                              "swift");
521

522
523
  /* How often to check for the stop file and dump restarts and exit the
   * application. */
524
  const int restart_stop_steps =
Peter W. Draper's avatar
Peter W. Draper committed
525
      parser_get_opt_param_int(params, "Restarts:stop_steps", 100);
526

527
528
529
530
531
532
533
534
535
536
537
  /* Get the maximal wall-clock time of this run */
  const float restart_max_hours_runtime =
      parser_get_opt_param_float(params, "Restarts:max_run_time", FLT_MAX);

  /* Do we want to resubmit when we hit the limit? */
  const int resubmit_after_max_hours =
      parser_get_opt_param_int(params, "Restarts:resubmit_on_exit", 0);

  /* What command should we run to resubmit at the end? */
  char resubmit_command[PARSER_MAX_LINE_SIZE];
  if (resubmit_after_max_hours)
538
539
    parser_get_param_string(params, "Restarts:resubmit_command",
                            resubmit_command);
540

541
  /* If restarting, look for the restart files. */
542
  if (restart) {
543

544
    /* Attempting a restart. */
545
546
    char **restart_files = NULL;
    int restart_nfiles = 0;
547
548

    if (myrank == 0) {
549
      message("Restarting SWIFT");
550

551
      /* Locate the restart files. */
Peter W. Draper's avatar
Peter W. Draper committed
552
553
      restart_files =
          restart_locate(restart_dir, restart_name, &restart_nfiles);
554
555
      if (restart_nfiles == 0)
        error("Failed to locate any restart files in %s", restart_dir);
556
557

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

      if (verbose > 0)
563
564
        for (int i = 0; i < restart_nfiles; i++)
          message("found restart file: %s", restart_files[i]);
565
566
567
    }

#ifdef WITH_MPI
568
    /* Distribute the restart files, need one for each rank. */
569
570
571
    if (myrank == 0) {

      for (int i = 1; i < nr_nodes; i++) {
572
573
        strcpy(restart_file, restart_files[i]);
        MPI_Send(restart_file, 200, MPI_BYTE, i, 0, MPI_COMM_WORLD);
574
575
576
      }

      /* Keep local file. */
577
      strcpy(restart_file, restart_files[0]);
578
579

      /* Finished with the list. */
580
      restart_locate_free(restart_nfiles, restart_files);
581
582

    } else {
583
      MPI_Recv(restart_file, 200, MPI_BYTE, 0, 0, MPI_COMM_WORLD,
584
585
               MPI_STATUS_IGNORE);
    }
586
    if (verbose > 1) message("local restart file = %s", restart_file);
587
#else
588

589
    /* Just one restart file. */
590
    strcpy(restart_file, restart_files[0]);
591
592
#endif

593
    /* Now read it. */
594
    restart_read(&e, restart_file);
595
596
597

    /* And initialize the engine with the space and policies. */
    if (myrank == 0) clocks_gettime(&tic);
Peter W. Draper's avatar
Peter W. Draper committed
598
599
    engine_config(1, &e, params, nr_nodes, myrank, nr_threads, with_aff,
                  talking, restart_file);
600
601
    if (myrank == 0) {
      clocks_gettime(&toc);
602
      message("engine_config took %.3f %s.", clocks_diff(&tic, &toc),
603
604
605
606
              clocks_getunit());
      fflush(stdout);
    }

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

611
  } else {
612

613
    /* Not restarting so look for the ICs. */
614
    /* Initialize unit system and constants */
615
    units_init_from_params(&us, params, "InternalUnitSystem");
616
    phys_const_init(&us, params, &prog_const);
617
    if (myrank == 0) {
618
619
620
621
622
623
624
      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
625

626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
    /* Read particles and space information from ICs */
    char ICfileName[200] = "";
    parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
    const int periodic =
        parser_get_param_int(params, "InitialConditions:periodic");
    const int replicate =
        parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
    clean_smoothing_length_values = parser_get_opt_param_int(
        params, "InitialConditions:cleanup_smoothing_lengths", 0);
    const int cleanup_h = parser_get_opt_param_int(
        params, "InitialConditions:cleanup_h_factors", 0);
    const int cleanup_sqrt_a = parser_get_opt_param_int(
        params, "InitialConditions:cleanup_velocity_factors", 0);
    const int generate_gas_in_ics = parser_get_opt_param_int(
        params, "InitialConditions:generate_gas_in_ics", 0);

    /* Some checks that we are not doing something stupid */
    if (generate_gas_in_ics && flag_entropy_ICs)
      error("Can't generate gas if the entropy flag is set in the ICs.");
    if (generate_gas_in_ics && !with_cosmology)
      error("Can't generate gas if the run is not cosmological.");

648
    /* Initialise the cosmology */
649
650
651
    if (with_cosmology)
      cosmology_init(params, &us, &prog_const, &cosmo);
    else
652
      cosmology_init_no_cosmo(&cosmo);
653
    if (myrank == 0 && with_cosmology) cosmology_print(&cosmo);
654

655
    /* Initialise the hydro properties */
656
657
    if (with_hydro)
      hydro_props_init(&hydro_properties, &prog_const, &us, params);
658
659
660
661
662
663
664
665
    else
      bzero(&hydro_properties, sizeof(struct hydro_props));

    /* Initialise the equation of state */
    if (with_hydro)
      eos_init(&eos, &prog_const, &us, params);
    else
      bzero(&eos, sizeof(struct eos_parameters));
666

Loic Hausammann's avatar
Loic Hausammann committed
667
668
    /* Initialise the stars properties */
    if (with_stars)
669
670
      stars_props_init(&stars_properties, &prog_const, &us, params,
                       &hydro_properties);
Loic Hausammann's avatar
Loic Hausammann committed
671
672
673
    else
      bzero(&stars_properties, sizeof(struct stars_props));

674
    /* Initialise the gravity properties */
675
    if (with_self_gravity)
676
677
      gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology,
                         periodic);
678
679
    else
      bzero(&gravity_properties, sizeof(struct gravity_props));
680

681
    /* Be verbose about what happens next */
682
    if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
683
684
    if (myrank == 0 && cleanup_h)
      message("Cleaning up h-factors (h=%f)", cosmo.h);
685
686
    if (myrank == 0 && cleanup_sqrt_a)
      message("Cleaning up a-factors from velocity (a=%f)", cosmo.a);
687
    fflush(stdout);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
688

689
690
691
692
    /* Get ready to read particles of all kinds */
    size_t Ngas = 0, Ngpart = 0, Nspart = 0;
    double dim[3] = {0., 0., 0.};
    if (myrank == 0) clocks_gettime(&tic);
693
#if defined(HAVE_HDF5)
694
695
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
696
    read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
697
                     &Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
698
                     (with_external_gravity || with_self_gravity), with_stars,
699
700
701
                     cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
                     nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
                     dry_run);
702
703
#else
    read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
704
                   &Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
Matthieu Schaller's avatar
Matthieu Schaller committed
705
                   (with_external_gravity || with_self_gravity), with_stars,
706
707
708
                   cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
                   nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
                   dry_run);
709
710
#endif
#else
711
    read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
712
                   &Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
713
                   (with_external_gravity || with_self_gravity), with_stars,
714
715
                   cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads,
                   dry_run);
716
#endif
717
#endif
718
719
720
721
722
723
    if (myrank == 0) {
      clocks_gettime(&toc);
      message("Reading initial conditions took %.3f %s.",
              clocks_diff(&tic, &toc), clocks_getunit());
      fflush(stdout);
    }
724

725
#ifdef SWIFT_DEBUG_CHECKS
726
    /* Check once and for all that we don't have unwanted links */
727
    if (!with_stars && !dry_run) {
728
      for (size_t k = 0; k < Ngpart; ++k)
729
        if (gparts[k].type == swift_type_stars) error("Linking problem");
730
    }
731
    if (!with_hydro && !dry_run) {
732
733
734
      for (size_t k = 0; k < Ngpart; ++k)
        if (gparts[k].type == swift_type_gas) error("Linking problem");
    }
735
736

    /* Check that the other links are correctly set */
737
738
    if (!dry_run)
      part_verify_links(parts, gparts, sparts, Ngas, Ngpart, Nspart, 1);
739
#endif
740

741
742
    /* Get the total number of particles across all nodes. */
    long long N_total[3] = {0, 0, 0};
743
#if defined(WITH_MPI)
744
745
746
    long long N_long[3] = {Ngas, Ngpart, Nspart};
    MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
                  MPI_COMM_WORLD);
747
#else
748
749
750
    N_total[0] = Ngas;
    N_total[1] = Ngpart;
    N_total[2] = Nspart;
751
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
752

753
754
    if (myrank == 0)
      message(
755
          "Read %lld gas particles, %lld stars particles and %lld gparts from "
Loic Hausammann's avatar
Loic Hausammann committed
756
          "the ICs.",
757
          N_total[0], N_total[2], N_total[1]);
758

759
760
761
    /* Verify that the fields to dump actually exist */
    if (myrank == 0) io_check_output_fields(params, N_total);

762
763
    /* Initialize the space with these data. */
    if (myrank == 0) clocks_gettime(&tic);
764
    space_init(&s, params, &cosmo, dim, parts, gparts, sparts, Ngas, Ngpart,
765
               Nspart, periodic, replicate, generate_gas_in_ics, with_hydro,
766
               with_self_gravity, talking, dry_run);
767
768
769
770
771
772
773
774

    if (myrank == 0) {
      clocks_gettime(&toc);
      message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
              clocks_getunit());
      fflush(stdout);
    }

775
776
777
    /* Initialise the long-range gravity mesh */
    if (with_self_gravity && periodic) {
#ifdef HAVE_FFTW
778
      pm_mesh_init(&mesh, &gravity_properties, s.dim, nr_threads);
779
780
781
782
783
784
785
786
787
#else
      /* Need the FFTW library if periodic and self gravity. */
      error(
          "No FFTW library found. Cannot compute periodic long-range forces.");
#endif
    } else {
      pm_mesh_init_no_mesh(&mesh, s.dim);
    }

788
789
    /* Check that the matter content matches the cosmology given in the
     * parameter file. */
790
    if (with_cosmology && with_self_gravity && !dry_run)
791
      space_check_cosmology(&s, &cosmo, myrank);
792

793
/* Also update the total counts (in case of changes due to replication) */
794
#if defined(WITH_MPI)
795
796
797
798
799
    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);
800
#else
801
802
803
    N_total[0] = s.nr_parts;
    N_total[1] = s.nr_gparts;
    N_total[2] = s.nr_sparts;
804
805
#endif

806
807
808
809
810
811
812
813
814
815
816
817
818
    /* 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);
    }
819

820
821
822
823
824
825
    /* 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);
    }
826

827
828
829
830
831
832
    /* 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]);
    }
833

834
835
836
837
    /* Initialise the external potential properties */
    if (with_external_gravity)
      potential_init(params, &prog_const, &us, &s, &potential);
    if (myrank == 0) potential_print(&potential);
838

839
840
    /* Initialise the cooling function properties */
    if (with_cooling) cooling_init(params, &us, &prog_const, &cooling_func);
841
    if (myrank == 0) cooling_print(&cooling_func);
842

843
844
845
846
    /* Initialise the chemistry */
    chemistry_init(params, &us, &prog_const, &chemistry);
    if (myrank == 0) chemistry_print(&chemistry);

847
848
    /* Initialise the feedback properties */
    if (with_sourceterms) sourceterms_init(params, &us, &sourceterms);
Peter W. Draper's avatar
Peter W. Draper committed
849
    if (with_sourceterms && myrank == 0) sourceterms_print(&sourceterms);
850
851
852
853
854
855
856
857
858
859
860
861
862
863

    /* 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;
Loic Hausammann's avatar
Loic Hausammann committed
864
    if (with_feedback) engine_policies |= engine_policy_feedback;
Matthieu Schaller's avatar
Matthieu Schaller committed
865
866
    if (with_structure_finding)
      engine_policies |= engine_policy_structure_finding;
867

868
869
870
    // MATTHIEU: Temporary star formation law
    // engine_policies |= engine_policy_star_formation;

871
872
    /* Initialize the engine with the space and policies. */
    if (myrank == 0) clocks_gettime(&tic);
873
874
    engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
                engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
875
876
                &hydro_properties, &gravity_properties, &stars_properties,
                &mesh, &potential, &cooling_func, &chemistry, &sourceterms);
Peter W. Draper's avatar
Peter W. Draper committed
877
878
    engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
                  talking, restart_file);
879

880
881
882
883
884
885
    if (myrank == 0) {
      clocks_gettime(&toc);
      message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
              clocks_getunit());
      fflush(stdout);
    }
886

887
888
889
890
    /* Get some info to the user. */
    if (myrank == 0) {
      long long N_DM = N_total[1] - N_total[2] - N_total[0];
      message(
891
          "Running on %lld gas particles, %lld stars particles and %lld DM "
892
893
894
          "particles (%lld gravity particles)",
          N_total[0], N_total[2], N_total[1] > 0 ? N_DM : 0, N_total[1]);
      message(
895
896
897
          "from t=%.3e until t=%.3e with %d ranks, %d threads / rank and %d "
          "task queues / rank (dt_min=%.3e, dt_max=%.3e)...",
          e.time_begin, e.time_end, nr_nodes, e.nr_threads, e.sched.nr_queues,
898
          e.dt_min, e.dt_max);
899
900
      fflush(stdout);
    }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
901
  }
902

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
903
904
  /* Time to say good-bye if this was not a serious run. */
  if (dry_run) {
905
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
906
907
    if ((res = MPI_Finalize()) != MPI_SUCCESS)
      error