main.c 40.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
void print_help_message(void) {
58

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
  printf("  %2s %14s %s\n", "-o", "{str}",
lhausamm's avatar
Format    
lhausamm committed
94
         "Generate a default output parameter file.");
95
96
97
  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.");
98
  printf("  %2s %14s %s\n", "-r", "", "Continue using restart files.");
99
100
  printf("  %2s %14s %s\n", "-s", "", "Run with hydrodynamics.");
  printf("  %2s %14s %s\n", "-S", "", "Run with stars.");
Loic Hausammann's avatar
Loic Hausammann committed
101
  printf("  %2s %14s %s\n", "-b", "", "Run with stars feedback.");
102
  printf("  %2s %14s %s\n", "-t", "{int}",
103
104
         "The number of threads to use on each MPI rank. Defaults to 1 if not "
         "specified.");
105
106
107
108
  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.");
109
  printf("  %2s %14s %s\n", "-x", "", "Run with structure finding.");
110
  printf("  %2s %14s %s\n", "-y", "{int}",
111
         "Time-step frequency at which task graphs are dumped.");
112
113
  printf("  %2s %14s %s\n", "-Y", "{int}",
         "Time-step frequency at which threadpool tasks are dumped.");
114
  printf("  %2s %14s %s\n", "-h", "", "Print this help message and exit.");
115
116
117
  printf(
      "\nSee the file parameter_example.yml for an example of "
      "parameter file.\n");
118
119
}

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

126
  struct clocks_time tic, toc;
127
  struct engine e;
128

129
130
  /* Structs used by the engine. Declare now to make sure these are always in
   * scope.  */
lhausamm's avatar
lhausamm committed
131
  struct chemistry_global_data chemistry;
132
  struct cooling_function_data cooling_func;
133
  struct cosmology cosmo;
134
  struct external_potential potential;
135
  struct pm_mesh mesh;
136
137
138
  struct gpart *gparts = NULL;
  struct gravity_props gravity_properties;
  struct hydro_props hydro_properties;
Loic Hausammann's avatar
Loic Hausammann committed
139
  struct stars_props stars_properties;
140
141
142
143
144
145
  struct part *parts = NULL;
  struct phys_const prog_const;
  struct sourceterms sourceterms;
  struct space s;
  struct spart *sparts = NULL;
  struct unit_system us;
146

147
  int nr_nodes = 1, myrank = 0;
148

149
#ifdef WITH_MPI
150
  /* Start by initializing MPI. */
151
  int res = 0, prov = 0;
152
153
  if ((res = MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &prov)) !=
      MPI_SUCCESS)
154
155
    error("Call to MPI_Init failed with error %i.", res);
  if (prov != MPI_THREAD_MULTIPLE)
156
    error(
157
158
        "MPI does not provide the level of threading"
        " required (MPI_THREAD_MULTIPLE).");
159
  if ((res = MPI_Comm_size(MPI_COMM_WORLD, &nr_nodes)) != MPI_SUCCESS)
160
161
162
    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);
163
164
165
166

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

167
168
169
  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);
170
171
  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
172
           nr_nodes);
173
174
175
176
  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.");
  }
177
  fflush(stdout);
178

179
#endif
180

181
182
  /* Welcome to SWIFT, you made the right choice */
  if (myrank == 0) greetings();
183

184
  int with_aff = 0;
185
  int dry_run = 0;
186
  int dump_tasks = 0;
187
  int dump_threadpool = 0;
188
  int nsteps = -2;
189
  int restart = 0;
190
191
  int with_cosmology = 0;
  int with_external_gravity = 0;
Tom Theuns's avatar
Tom Theuns committed
192
  int with_sourceterms = 0;
193
  int with_cooling = 0;
194
195
  int with_self_gravity = 0;
  int with_hydro = 0;
196
  int with_stars = 0;
Loic Hausammann's avatar
Loic Hausammann committed
197
  int with_feedback = 0;
198
  int with_fp_exceptions = 0;
199
  int with_drift_all = 0;
200
  int with_mpole_reconstruction = 0;
201
  int with_structure_finding = 0;
202
  int verbose = 0;
203
  int nr_threads = 1;
204
  int with_verbose_timers = 0;
205
  int nparams = 0;
206
  char output_parameters_filename[200] = "";
207
  char *cmdparams[PARSER_MAX_NO_OF_PARAMS];
208
  char paramFileName[200] = "";
209
  char restart_file[200] = "";
210
211
212
  unsigned long long cpufreq = 0;

  /* Parse the parameters */
213
  int c;
Loic Hausammann's avatar
Loic Hausammann committed
214
  while ((c = getopt(argc, argv, "abcCdDef:FgGhMn:o:P:rsSt:Tv:xy:Y:")) != -1)
215
    switch (c) {
216
      case 'a':
217
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA)
218
        with_aff = 1;
219
220
221
#else
        error("Need NUMA support for thread affinity");
#endif
222
        break;
Loic Hausammann's avatar
Loic Hausammann committed
223
224
225
      case 'b':
        with_feedback = 1;
        break;
226
      case 'c':
227
        with_cosmology = 1;
228
        break;
229
230
231
      case 'C':
        with_cooling = 1;
        break;
232
      case 'd':
233
        dry_run = 1;
234
        break;
235
236
237
      case 'D':
        with_drift_all = 1;
        break;
238
      case 'e':
239
#ifdef HAVE_FE_ENABLE_EXCEPT
240
        with_fp_exceptions = 1;
241
242
243
#else
        error("Need support for floating point exception on this platform");
#endif
244
        break;
245
      case 'f':
246
247
248
        if (sscanf(optarg, "%llu", &cpufreq) != 1) {
          if (myrank == 0) printf("Error parsing CPU frequency (-f).\n");
          if (myrank == 0) print_help_message();
249
          return 1;
250
        }
251
        break;
252
253
254
      case 'F':
        with_sourceterms = 1;
        break;
255
      case 'g':
256
        with_external_gravity = 1;
257
        break;
258
259
      case 'G':
        with_self_gravity = 1;
Tom Theuns's avatar
Tom Theuns committed
260
        break;
261
      case 'h':
262
        if (myrank == 0) print_help_message();
263
        return 0;
264
265
266
      case 'M':
        with_mpole_reconstruction = 1;
        break;
267
268
269
270
271
272
273
      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;
274
      case 'o':
lhausamm's avatar
Format    
lhausamm committed
275
276
277
278
279
280
281
282
        if (sscanf(optarg, "%s", output_parameters_filename) != 1) {
          if (myrank == 0) {
            printf("Error parsing output fields filename");
            print_help_message();
          }
          return 1;
        }
        break;
283
284
285
286
      case 'P':
        cmdparams[nparams] = optarg;
        nparams++;
        break;
287
288
289
      case 'r':
        restart = 1;
        break;
290
      case 's':
291
        with_hydro = 1;
292
        break;
293
294
295
      case 'S':
        with_stars = 1;
        break;
296
297
298
299
300
301
302
303
      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;
304
305
306
      case 'T':
        with_verbose_timers = 1;
        break;
307
      case 'v':
308
309
310
        if (sscanf(optarg, "%d", &verbose) != 1) {
          if (myrank == 0) printf("Error parsing verbosity level (-v).\n");
          if (myrank == 0) print_help_message();
311
          return 1;
312
        }
313
        break;
314
      case 'x':
315
#ifdef HAVE_VELOCIRAPTOR
316
        with_structure_finding = 1;
317
#else
Matthieu Schaller's avatar
Matthieu Schaller committed
318
319
320
        error(
            "Error: (-x) needs to have the code compiled with VELOCIraptor "
            "linked in.");
321
#endif
322
        break;
323
      case 'y':
324
325
326
        if (sscanf(optarg, "%d", &dump_tasks) != 1) {
          if (myrank == 0) printf("Error parsing dump_tasks (-y). \n");
          if (myrank == 0) print_help_message();
327
          return 1;
328
        }
329
330
331
332
333
334
#ifndef SWIFT_DEBUG_TASKS
        if (dump_tasks) {
          error(
              "Task dumping is only possible if SWIFT was configured with the "
              "--enable-task-debugging option.");
        }
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
#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.");
        }
350
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
351
352
353
354
355
        break;
      case '?':
        if (myrank == 0) print_help_message();
        return 1;
        break;
356
    }
357
358

  /* Write output parameter file */
359
  if (myrank == 0 && strcmp(output_parameters_filename, "") != 0) {
360
    io_write_output_field_parameter(output_parameters_filename);
361
    printf("End of run.\n");
362
363
364
365
    return 0;
  }

  /* check inputs */
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
366
367
368
369
370
  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");
371
372
    if (myrank == 0) print_help_message();
    return 1;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
373
374
375
376
377
378
379
380
381
382
383
  } 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;
  }
384
385
386
  if (with_stars && !with_external_gravity && !with_self_gravity) {
    if (myrank == 0)
      printf(
387
          "Error: Cannot process stars without gravity, -g or -G must be "
388
389
390
391
          "chosen.\n");
    if (myrank == 0) print_help_message();
    return 1;
  }
392

Loic Hausammann's avatar
Loic Hausammann committed
393
394
395
396
397
398
399
400
401
  if (!with_stars && with_feedback) {
    if (myrank == 0)
      printf(
          "Error: Cannot process feedback without stars, -S must be "
          "chosen.\n");
    if (myrank == 0) print_help_message();
    return 1;
  }

402
403
/* Let's pin the main thread, now we know if affinity will be used. */
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
404
405
  if (with_aff &&
      ((ENGINE_POLICY)&engine_policy_setaffinity) == engine_policy_setaffinity)
406
407
408
    engine_pin();
#endif

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

412
413
414
  /* How vocal are we ? */
  const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
419
420
421
422
423
  /* Report CPU frequency.*/
  cpufreq = clocks_get_cpufreq();
  if (myrank == 0) {
    message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
  }
424

Matthieu Schaller's avatar
Matthieu Schaller committed
425
/* Report host name(s). */
426
#ifdef WITH_MPI
427
  if (talking) {
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
428
429
    message("Rank %d running on: %s", myrank, hostname());
  }
430
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
431
  message("Running on: %s", hostname());
432
433
#endif

434
435
/* Do we have debugging checks ? */
#ifdef SWIFT_DEBUG_CHECKS
436
437
  if (myrank == 0)
    message("WARNING: Debugging checks activated. Code will be slower !");
438
439
#endif

440
441
442
443
444
445
446
/* Do we have debugging checks ? */
#ifdef SWIFT_USE_NAIVE_INTERACTIONS
  if (myrank == 0)
    message(
        "WARNING: Naive cell interactions activated. Code will be slower !");
#endif

447
448
449
450
451
452
453
454
455
/* 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
456
457
  /* Do we choke on FP-exceptions ? */
  if (with_fp_exceptions) {
458
#ifdef HAVE_FE_ENABLE_EXCEPT
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
459
    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
460
#endif
461
462
    if (myrank == 0)
      message("WARNING: Floating point exceptions will be reported.");
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
463
  }
464

465
/* Do we have slow barriers? */
466
#ifndef HAVE_PTHREAD_BARRIERS
467
  if (myrank == 0)
468
    message("WARNING: Non-optimal thread barriers are being used.");
469
470
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
471
472
  /* How large are the parts? */
  if (myrank == 0) {
473
474
475
476
477
478
479
480
    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
481
  }
482

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
483
  /* Read the parameter file */
Matthieu Schaller's avatar
Matthieu Schaller committed
484
485
  struct swift_params *params =
      (struct swift_params *)malloc(sizeof(struct swift_params));
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
486
487
488
489
  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);
490
491

    /* Handle any command-line overrides. */
492
493
494
495
    if (nparams > 0) {
      message(
          "Overwriting values read from the YAML file with command-line "
          "values.");
496
      for (int k = 0; k < nparams; k++) parser_set_param(params, cmdparams[k]);
497
    }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
498
  }
499
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
500
501
  /* Broadcast the parameter file */
  MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
502
#endif
503

504
505
  /* Temporary early aborts for modes not supported over MPI. */
#ifdef WITH_MPI
Matthieu Schaller's avatar
Matthieu Schaller committed
506
507
  if (with_mpole_reconstruction && nr_nodes > 1)
    error("Cannot reconstruct m-poles every step over MPI (yet).");
508
509
#endif

510
511
512
513
#ifdef WITH_MPI
  if (with_feedback) error("Can't run with feedback over MPI (yet).");
#endif

514
#if defined(WITH_MPI) && defined(HAVE_VELOCIRAPTOR)
Matthieu Schaller's avatar
Matthieu Schaller committed
515
516
  if (with_structure_finding && nr_nodes > 1)
    error("VEOCIraptor not yet enabled over MPI.");
517
518
#endif

519
  /* Check that we can write the snapshots by testing if the output
520
   * directory exists and is searchable and writable. */
521
522
523
  char basename[PARSER_MAX_LINE_SIZE];
  parser_get_param_string(params, "Snapshots:basename", basename);
  const char *dirp = dirname(basename);
524
525
  if (access(dirp, W_OK | X_OK) != 0) {
    error("Cannot write snapshots in directory %s (%s)", dirp, strerror(errno));
526
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
527
528
529

  /* Check that we can write the structure finding catalogues by testing if the
   * output
530
   * directory exists and is searchable and writable. */
Matthieu Schaller's avatar
Matthieu Schaller committed
531
  if (with_structure_finding) {
532
    char stfbasename[PARSER_MAX_LINE_SIZE];
533
    parser_get_param_string(params, "StructureFinding:basename", stfbasename);
534
535
    const char *stfdirp = dirname(stfbasename);
    if (access(stfdirp, W_OK | X_OK) != 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
536
537
      error("Cannot write stf catalogues in directory %s (%s)", stfdirp,
            strerror(errno));
538
539
    }
  }
540

541
  /* Prepare the domain decomposition scheme */
542
  struct repartition reparttype;
543
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
544
545
546
547
548
  struct partition initial_partition;
  partition_init(&initial_partition, &reparttype, params, nr_nodes);

  /* Let's report what we did */
  if (myrank == 0) {
549
550
#if defined(HAVE_PARMETIS)
    if (reparttype.usemetis)
Matthieu Schaller's avatar
Matthieu Schaller committed
551
      message("Using METIS serial partitioning:");
552
    else
Matthieu Schaller's avatar
Matthieu Schaller committed
553
      message("Using ParMETIS partitioning:");
554
555
556
557
#else
    message("Using METIS serial partitioning:");
#endif
    message("  initial partitioning: %s",
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
558
559
            initial_partition_name[initial_partition.type]);
    if (initial_partition.type == INITPART_GRID)
560
      message("    grid set to [ %i %i %i ].", initial_partition.grid[0],
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
561
              initial_partition.grid[1], initial_partition.grid[2]);
562
    message("  repartitioning: %s", repartition_name[reparttype.type]);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
563
  }
564
#endif
565

566
  /* Common variables for restart and IC sections. */
567
  int clean_smoothing_length_values = 0;
568
569
570
  int flag_entropy_ICs = 0;

  /* Work out where we will read and write restart files. */
571
  char restart_dir[PARSER_MAX_LINE_SIZE];
Peter W. Draper's avatar
Peter W. Draper committed
572
573
  parser_get_opt_param_string(params, "Restarts:subdir", restart_dir,
                              "restart");
574
575

  /* The directory must exist. */
576
  if (myrank == 0) {
577
    if (access(restart_dir, W_OK | X_OK) != 0) {
578
      if (restart) {
579
        error("Cannot restart as no restart subdirectory: %s (%s)", restart_dir,
580
              strerror(errno));
581
      } else {
Peter W. Draper's avatar
Peter W. Draper committed
582
583
584
        if (mkdir(restart_dir, 0777) != 0)
          error("Failed to create restart directory: %s (%s)", restart_dir,
                strerror(errno));
585
      }
586
    }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
587
588
  }

589
  /* Basename for any restart files. */
590
591
  char restart_name[PARSER_MAX_LINE_SIZE];
  parser_get_opt_param_string(params, "Restarts:basename", restart_name,
592
                              "swift");
593

594
595
  /* How often to check for the stop file and dump restarts and exit the
   * application. */
596
  const int restart_stop_steps =
Peter W. Draper's avatar
Peter W. Draper committed
597
      parser_get_opt_param_int(params, "Restarts:stop_steps", 100);
598

599
600
601
602
603
604
605
606
607
608
609
  /* 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)
610
611
    parser_get_param_string(params, "Restarts:resubmit_command",
                            resubmit_command);
612

613
  /* If restarting, look for the restart files. */
614
  if (restart) {
615

616
    /* Attempting a restart. */
617
618
    char **restart_files = NULL;
    int restart_nfiles = 0;
619
620

    if (myrank == 0) {
621
      message("Restarting SWIFT");
622

623
      /* Locate the restart files. */
Peter W. Draper's avatar
Peter W. Draper committed
624
625
      restart_files =
          restart_locate(restart_dir, restart_name, &restart_nfiles);
626
627
      if (restart_nfiles == 0)
        error("Failed to locate any restart files in %s", restart_dir);
628
629

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

      if (verbose > 0)
635
636
        for (int i = 0; i < restart_nfiles; i++)
          message("found restart file: %s", restart_files[i]);
637
638
639
    }

#ifdef WITH_MPI
640
    /* Distribute the restart files, need one for each rank. */
641
642
643
    if (myrank == 0) {

      for (int i = 1; i < nr_nodes; i++) {
644
645
        strcpy(restart_file, restart_files[i]);
        MPI_Send(restart_file, 200, MPI_BYTE, i, 0, MPI_COMM_WORLD);
646
647
648
      }

      /* Keep local file. */
649
      strcpy(restart_file, restart_files[0]);
650
651

      /* Finished with the list. */
652
      restart_locate_free(restart_nfiles, restart_files);
653
654

    } else {
655
      MPI_Recv(restart_file, 200, MPI_BYTE, 0, 0, MPI_COMM_WORLD,
656
657
               MPI_STATUS_IGNORE);
    }
658
    if (verbose > 1) message("local restart file = %s", restart_file);
659
#else
660

661
    /* Just one restart file. */
662
    strcpy(restart_file, restart_files[0]);
663
664
#endif

665
    /* Now read it. */
666
    restart_read(&e, restart_file);
667
668
669

    /* And initialize the engine with the space and policies. */
    if (myrank == 0) clocks_gettime(&tic);
Peter W. Draper's avatar
Peter W. Draper committed
670
671
    engine_config(1, &e, params, nr_nodes, myrank, nr_threads, with_aff,
                  talking, restart_file);
672
673
    if (myrank == 0) {
      clocks_gettime(&toc);
674
      message("engine_config took %.3f %s.", clocks_diff(&tic, &toc),
675
676
677
678
              clocks_getunit());
      fflush(stdout);
    }

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

683
  } else {
684

685
    /* Not restarting so look for the ICs. */
686
    /* Initialize unit system and constants */
687
    units_init_from_params(&us, params, "InternalUnitSystem");
688
    phys_const_init(&us, params, &prog_const);
689
    if (myrank == 0) {
690
691
692
693
694
695
696
      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
697

698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
    /* 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.");

720
    /* Initialise the cosmology */
721
722
723
    if (with_cosmology)
      cosmology_init(params, &us, &prog_const, &cosmo);
    else
724
      cosmology_init_no_cosmo(&cosmo);
725
    if (myrank == 0 && with_cosmology) cosmology_print(&cosmo);
726

727
    /* Initialise the hydro properties */
728
729
    if (with_hydro)
      hydro_props_init(&hydro_properties, &prog_const, &us, params);
730
731
732
733
734
735
736
737
    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));
738

Loic Hausammann's avatar
Loic Hausammann committed
739
740
    /* Initialise the stars properties */
    if (with_stars)
741
742
      stars_props_init(&stars_properties, &prog_const, &us, params,
                       &hydro_properties);
Loic Hausammann's avatar
Loic Hausammann committed
743
744
745
    else
      bzero(&stars_properties, sizeof(struct stars_props));

746
    /* Initialise the gravity properties */
747
    if (with_self_gravity)
748
749
      gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology,
                         periodic);
750
751
    else
      bzero(&gravity_properties, sizeof(struct gravity_props));
752

753
    /* Be verbose about what happens next */
754
    if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
755
756
    if (myrank == 0 && cleanup_h)
      message("Cleaning up h-factors (h=%f)", cosmo.h);
757
758
    if (myrank == 0 && cleanup_sqrt_a)
      message("Cleaning up a-factors from velocity (a=%f)", cosmo.a);
759
    fflush(stdout);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
760

761
762
763
764
    /* 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);
765
#if defined(HAVE_HDF5)
766
767
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
768
    read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
769
                     &Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
770
                     (with_external_gravity || with_self_gravity), with_stars,
771
772
773
                     cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
                     nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
                     dry_run);
774
775
#else
    read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
776
                   &Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
Matthieu Schaller's avatar
Matthieu Schaller committed
777
                   (with_external_gravity || with_self_gravity), with_stars,
778
779
780
                   cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
                   nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
                   dry_run);
781
782
#endif
#else
783
    read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
784
                   &Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
785
                   (with_external_gravity || with_self_gravity), with_stars,
786
787
                   cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads,
                   dry_run);
788
#endif
789
#endif
790
791
792
793
794
795
    if (myrank == 0) {
      clocks_gettime(&toc);
      message("Reading initial conditions took %.3f %s.",
              clocks_diff(&tic, &toc), clocks_getunit());
      fflush(stdout);
    }
796

797
#ifdef SWIFT_DEBUG_CHECKS
798
    /* Check once and for all that we don't have unwanted links */
799
    if (!with_stars && !dry_run) {
800
      for (size_t k = 0; k < Ngpart; ++k)
801
        if (gparts[k].type == swift_type_stars) error("Linking problem");
802
    }
803
    if (!with_hydro && !dry_run) {
804
805
806
      for (size_t k = 0; k < Ngpart; ++k)
        if (gparts[k].type == swift_type_gas) error("Linking problem");
    }
807
808

    /* Check that the other links are correctly set */
809
810
    if (!dry_run)
      part_verify_links(parts, gparts, sparts, Ngas, Ngpart, Nspart, 1);
811
#endif
812

813
814
    /* Get the total number of particles across all nodes. */
    long long N_total[3] = {0, 0, 0};
815
#if defined(WITH_MPI)
816
817
818
    long long N_long[3] = {Ngas, Ngpart, Nspart};
    MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
                  MPI_COMM_WORLD);
819
#else
820
821
822
    N_total[0] = Ngas;
    N_total[1] = Ngpart;
    N_total[2] = Nspart;
823
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
824

825
826
    if (myrank == 0)
      message(
827
          "Read %lld gas particles, %lld stars particles and %lld gparts from "
Loic Hausammann's avatar
Loic Hausammann committed
828
          "the ICs.",
829
          N_total[0], N_total[2], N_total[1]);
830

831
832
833
    /* Verify that the fields to dump actually exist */
    if (myrank == 0) io_check_output_fields(params, N_total);

834