main.c 23.7 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 <fenv.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
30
31
32
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
33
#include <unistd.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
34

35
36
/* MPI headers. */
#ifdef WITH_MPI
37
#include <mpi.h>
38
39
#endif

Pedro Gonnet's avatar
Pedro Gonnet committed
40
/* Local headers. */
41
#include "swift.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
42

43
44
/* Engine policy flags. */
#ifndef ENGINE_POLICY
45
#define ENGINE_POLICY engine_policy_none
46
47
#endif

James Willis's avatar
James Willis committed
48
49
50
/* Global profiler. */
struct profiler prof;

51
52
53
/**
 * @brief Help messages for the command line parameters.
 */
54
55
void print_help_message() {

Matthieu Schaller's avatar
Matthieu Schaller committed
56
57
  printf("\nUsage: swift [OPTION]... PARAMFILE\n");
  printf("       swift_mpi [OPTION]... PARAMFILE\n");
58

59
  printf("Valid options are:\n");
60
  printf("  %2s %8s %s\n", "-a", "", "Pin runners using processor affinity");
61
  printf("  %2s %8s %s\n", "-c", "", "Run with cosmological time integration");
62
  printf("  %2s %8s %s\n", "-C", "", "Run with cooling");
63
64
65
66
67
68
69
70
71
  printf(
      "  %2s %8s %s\n", "-d", "",
      "Dry run. Read the parameter file, allocate memory but does not read ");
  printf(
      "  %2s %8s %s\n", "", "",
      "the particles from ICs and exit before the start of time integration.");
  printf("  %2s %8s %s\n", "", "",
         "Allows user to check validy of parameter and IC files as well as "
         "memory limits.");
72
73
  printf("  %2s %8s %s\n", "-D", "",
         "Always drift all particles even the ones far from active particles.");
74
75
  printf("  %2s %8s %s\n", "-e", "",
         "Enable floating-point exceptions (debugging mode)");
76
77
  printf("  %2s %8s %s\n", "-f", "{int}",
         "Overwrite the CPU frequency (Hz) to be used for time measurements");
78
79
  printf("  %2s %8s %s\n", "-g", "",
         "Run with an external gravitational potential");
Tom Theuns's avatar
Tom Theuns committed
80
  printf("  %2s %8s %s\n", "-F", "", "Run with feedback ");
81
  printf("  %2s %8s %s\n", "-G", "", "Run with self-gravity");
Peter W. Draper's avatar
Peter W. Draper committed
82
  printf("  %2s %8s %s\n", "-n", "{int}",
83
84
         "Execute a fixed number of time steps. When unset use the time_end "
         "parameter to stop.");
85
  printf("  %2s %8s %s\n", "-s", "", "Run with SPH");
86
  printf("  %2s %8s %s\n", "-S", "", "Run with stars");
87
88
89
  printf("  %2s %8s %s\n", "-t", "{int}",
         "The number of threads to use on each MPI rank. Defaults to 1 if not "
         "specified.");
Matthieu Schaller's avatar
Matthieu Schaller committed
90
91
  printf("  %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity");
  printf("  %2s %8s %s\n", "", "", "1: MPI-rank 0 writes ");
92
  printf("  %2s %8s %s\n", "", "", "2: All MPI-ranks write");
93
  printf("  %2s %8s %s\n", "-y", "{int}",
94
         "Time-step frequency at which task graphs are dumped");
95
  printf("  %2s %8s %s\n", "-h", "", "Print this help message and exit");
96
97
98
  printf(
      "\nSee the file parameter_example.yml for an example of "
      "parameter file.\n");
99
100
}

Pedro Gonnet's avatar
Pedro Gonnet committed
101
102
103
104
/**
 * @brief Main routine that loads a few particles and generates some output.
 *
 */
105
106
int main(int argc, char *argv[]) {

107
  struct clocks_time tic, toc;
108

109
  int nr_nodes = 1, myrank = 0;
110

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

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

129
130
131
  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);
132
  if (myrank == 0) message("MPI is up and running with %i node(s).", nr_nodes);
133
134
135
136
  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.");
  }
137
  fflush(stdout);
138

139
#endif
140

141
/* Let's pin the main thread */
142
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
143
  if (((ENGINE_POLICY)&engine_policy_setaffinity) == engine_policy_setaffinity)
144
    engine_pin();
Angus Lepper's avatar
Angus Lepper committed
145
146
#endif

147
148
  /* Welcome to SWIFT, you made the right choice */
  if (myrank == 0) greetings();
149

150
  int with_aff = 0;
151
  int dry_run = 0;
152
  int dump_tasks = 0;
153
  int nsteps = -2;
154
155
  int with_cosmology = 0;
  int with_external_gravity = 0;
Tom Theuns's avatar
Tom Theuns committed
156
  int with_sourceterms = 0;
157
  int with_cooling = 0;
158
159
  int with_self_gravity = 0;
  int with_hydro = 0;
160
  int with_stars = 0;
161
  int with_fp_exceptions = 0;
162
  int with_drift_all = 0;
163
  int verbose = 0;
164
  int nr_threads = 1;
165
166
167
168
  char paramFileName[200] = "";
  unsigned long long cpufreq = 0;

  /* Parse the parameters */
169
  int c;
170
  while ((c = getopt(argc, argv, "acCdDef:FgGhn:sSt:v:y:")) != -1) switch (c) {
171
      case 'a':
172
        with_aff = 1;
173
        break;
174
      case 'c':
175
        with_cosmology = 1;
176
        break;
177
178
179
      case 'C':
        with_cooling = 1;
        break;
180
      case 'd':
181
        dry_run = 1;
182
        break;
183
184
185
      case 'D':
        with_drift_all = 1;
        break;
186
      case 'e':
187
        with_fp_exceptions = 1;
188
        break;
189
      case 'f':
190
191
192
        if (sscanf(optarg, "%llu", &cpufreq) != 1) {
          if (myrank == 0) printf("Error parsing CPU frequency (-f).\n");
          if (myrank == 0) print_help_message();
193
          return 1;
194
        }
195
        break;
196
197
198
      case 'F':
        with_sourceterms = 1;
        break;
199
      case 'g':
200
        with_external_gravity = 1;
201
        break;
202
203
      case 'G':
        with_self_gravity = 1;
Tom Theuns's avatar
Tom Theuns committed
204
        break;
205
      case 'h':
206
        if (myrank == 0) print_help_message();
207
        return 0;
208
209
210
211
212
213
214
      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;
215
      case 's':
216
        with_hydro = 1;
217
        break;
218
219
220
      case 'S':
        with_stars = 1;
        break;
221
222
223
224
225
226
227
228
      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;
229
      case 'v':
230
231
232
        if (sscanf(optarg, "%d", &verbose) != 1) {
          if (myrank == 0) printf("Error parsing verbosity level (-v).\n");
          if (myrank == 0) print_help_message();
233
          return 1;
234
        }
235
        break;
236
      case 'y':
237
238
239
        if (sscanf(optarg, "%d", &dump_tasks) != 1) {
          if (myrank == 0) printf("Error parsing dump_tasks (-y). \n");
          if (myrank == 0) print_help_message();
240
          return 1;
241
        }
242
243
244
245
246
247
248
#ifndef SWIFT_DEBUG_TASKS
        if (dump_tasks) {
          error(
              "Task dumping is only possible if SWIFT was configured with the "
              "--enable-task-debugging option.");
        }
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
249
250
251
252
253
        break;
      case '?':
        if (myrank == 0) print_help_message();
        return 1;
        break;
254
    }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
255
256
257
258
259
  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");
260
261
    if (myrank == 0) print_help_message();
    return 1;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
262
263
264
265
266
267
268
269
270
271
272
  } 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;
  }
273

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

277
278
279
  /* How vocal are we ? */
  const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
284
285
286
287
288
  /* Report CPU frequency.*/
  cpufreq = clocks_get_cpufreq();
  if (myrank == 0) {
    message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
  }
289

Matthieu Schaller's avatar
Matthieu Schaller committed
290
/* Report host name(s). */
291
#ifdef WITH_MPI
292
  if (talking) {
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
293
294
    message("Rank %d running on: %s", myrank, hostname());
  }
295
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
296
  message("Running on: %s", hostname());
297
298
#endif

299
300
/* Do we have debugging checks ? */
#ifdef SWIFT_DEBUG_CHECKS
301
302
  if (myrank == 0)
    message("WARNING: Debugging checks activated. Code will be slower !");
303
304
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
305
306
307
  /* Do we choke on FP-exceptions ? */
  if (with_fp_exceptions) {
    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
308
309
    if (myrank == 0)
      message("WARNING: Floating point exceptions will be reported.");
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
310
  }
311

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
312
313
314
315
  /* How large are the parts? */
  if (myrank == 0) {
    message("sizeof(struct part)  is %4zi bytes.", sizeof(struct part));
    message("sizeof(struct xpart) is %4zi bytes.", sizeof(struct xpart));
316
    message("sizeof(struct spart) is %4zi bytes.", sizeof(struct spart));
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
317
318
319
320
    message("sizeof(struct gpart) is %4zi bytes.", sizeof(struct gpart));
    message("sizeof(struct task)  is %4zi bytes.", sizeof(struct task));
    message("sizeof(struct cell)  is %4zi bytes.", sizeof(struct cell));
  }
321

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
322
323
324
325
326
327
328
329
330
  /* 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);
    // parser_print_params(&params);
    parser_write_params_to_file(params, "used_parameters.yml");
  }
331
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
332
333
  /* Broadcast the parameter file */
  MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
334
#endif
335

336
  /* Prepare the domain decomposition scheme */
337
  enum repartition_type reparttype = REPART_NONE;
338
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
339
340
341
342
343
344
345
346
347
348
349
350
  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]);
    message("Using %s repartitioning", repartition_name[reparttype]);
  }
351
#endif
352

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
  /* Initialize unit system and constants */
  struct UnitSystem 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);
  }

  /* Initialise the hydro properties */
  struct hydro_props hydro_properties;
Matthieu Schaller's avatar
Matthieu Schaller committed
369
  if (with_hydro) hydro_props_init(&hydro_properties, params);
370

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
371
372
373
374
375
376
  /* Read particles and space information from (GADGET) ICs */
  char ICfileName[200] = "";
  parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
  if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
  fflush(stdout);

377
  /* Get ready to read particles of all kinds */
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
378
379
  struct part *parts = NULL;
  struct gpart *gparts = NULL;
380
381
  struct spart *sparts = NULL;
  size_t Ngas = 0, Ngpart = 0, Nspart = 0;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
382
383
384
385
  double dim[3] = {0., 0., 0.};
  int periodic = 0;
  int flag_entropy_ICs = 0;
  if (myrank == 0) clocks_gettime(&tic);
386
387
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
Matthieu Schaller's avatar
Matthieu Schaller committed
388
389
390
391
  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, dry_run);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
392
#else
393
  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
Matthieu Schaller's avatar
Typos    
Matthieu Schaller committed
394
                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
395
396
                 (with_external_gravity || with_self_gravity), with_stars,
                 myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
397
398
#endif
#else
399
  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
400
401
402
                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                 (with_external_gravity || with_self_gravity), with_stars,
                 dry_run);
403
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
404
405
406
407
408
409
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("Reading initial conditions took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
410

411
412
#ifdef SWIFT_DEBUG_CHECKS
  /* Check once and for all that we don't have unwanted links */
413
  if (!with_stars) {
414
    for (size_t k = 0; k < Ngpart; ++k)
415
416
      if (gparts[k].type == swift_type_star) error("Linking problem");
  }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
417
418
  if (!with_hydro) {
    for (size_t k = 0; k < Ngpart; ++k)
419
      if (gparts[k].type == swift_type_gas) error("Linking problem");
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
420
  }
421
#endif
422

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
423
  /* Get the total number of particles across all nodes. */
424
  long long N_total[3] = {0, 0, 0};
425
#if defined(WITH_MPI)
426
  long long N_long[3] = {Ngas, Ngpart, Nspart};
Matthieu Schaller's avatar
Matthieu Schaller committed
427
428
  MPI_Reduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM, 0,
             MPI_COMM_WORLD);
429
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
430
431
  N_total[0] = Ngas;
  N_total[1] = Ngpart;
432
  N_total[2] = Nspart;
433
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
434
  if (myrank == 0)
435
436
437
438
    message(
        "Read %lld gas particles, %lld star particles and %lld gparts from the "
        "ICs.",
        N_total[0], N_total[2], N_total[1]);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
439
440
441
442

  /* Initialize the space with these data. */
  if (myrank == 0) clocks_gettime(&tic);
  struct space s;
443
444
  space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
             periodic, with_self_gravity, talking, dry_run);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
445
446
447
448
449
450
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
451

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
452
453
454
455
456
457
458
459
460
461
462
463
  /* 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("maximum depth is %d.", s.maxdepth);
    fflush(stdout);
  }
464

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
465
466
467
468
469
470
  /* 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);
  }
471

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
472
473
474
475
476
477
  /* 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]);
  }
478

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
479
480
481
  /* Initialise the external potential properties */
  struct external_potential potential;
  if (with_external_gravity)
482
    potential_init(params, &prog_const, &us, &s, &potential);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
  if (with_external_gravity && myrank == 0) potential_print(&potential);

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

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

  /* Construct the engine policy */
  int engine_policies = ENGINE_POLICY | engine_policy_steal;
  if (with_drift_all) engine_policies |= engine_policy_drift_all;
  if (with_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;
504
  if (with_stars) engine_policies |= engine_policy_stars;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
505
506
507
508
509

  /* Initialize the engine with the space and policies. */
  if (myrank == 0) clocks_gettime(&tic);
  struct engine e;
  engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
510
511
              engine_policies, talking, reparttype, &us, &prog_const,
              &hydro_properties, &potential, &cooling_func, &sourceterms);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
512
513
514
515
516
517
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
518

519
520
/* Init the runner history. */
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
521
  for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
522
523
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
524
525
  /* Get some info to the user. */
  if (myrank == 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
526
    long long N_DM = N_total[1] - N_total[2] - N_total[0];
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
527
    message(
Matthieu Schaller's avatar
Matthieu Schaller committed
528
529
530
531
532
533
534
535
        "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);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
536
537
    fflush(stdout);
  }
538

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
539
540
  /* Time to say good-bye if this was not a serious run. */
  if (dry_run) {
541
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
542
543
    if ((res = MPI_Finalize()) != MPI_SUCCESS)
      error("call to MPI_Finalize failed with error %i.", res);
544
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
545
546
547
548
549
550
    if (myrank == 0)
      message("Time integration ready to start. End of dry-run.");
    engine_clean(&e);
    free(params);
    return 0;
  }
551
552

#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
553
554
555
  /* Split the space. */
  engine_split(&e, &initial_partition);
  engine_redistribute(&e);
556
557
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
558
559
  /* Initialise the particles */
  engine_init_particles(&e, flag_entropy_ICs);
560

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
561
562
  /* Write the state of the system before starting time integration. */
  engine_dump_snapshot(&e);
563

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
564
565
  /* Legend */
  if (myrank == 0)
566
567
568
    printf("# %6s %14s %14s %10s %10s %10s %16s [%s]\n", "Step", "Time",
           "Time-step", "Updates", "g-Updates", "s-Updates", "Wall-clock time",
           clocks_getunit());
569

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
570
571
  /* Main simulation loop */
  for (int j = 0; !engine_is_done(&e) && e.step != nsteps; j++) {
572

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
573
574
    /* Reset timers */
    timers_reset(timers_mask_all);
575

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
576
577
    /* Take a step. */
    engine_step(&e);
578

579
#ifdef SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
580
581
    /* Dump the task data using the given frequency. */
    if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
582
583
#ifdef WITH_MPI

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
584
585
586
587
588
589
      /* 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");
590
        fclose(file_thread);
591
      }
592
      MPI_Barrier(MPI_COMM_WORLD);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633

      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");

          fprintf(file_thread, " %03i 0 0 0 0 %lli %lli 0 0 0 0 %lli\n", myrank,
                  e.tic_step, e.toc_step, cpufreq);
          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(
                  file_thread, " %03i %i %i %i %i %lli %lli %i %i %i %i %i\n",
                  myrank, e.sched.tasks[l].rid, e.sched.tasks[l].type,
                  e.sched.tasks[l].subtype, (e.sched.tasks[l].cj == NULL),
                  e.sched.tasks[l].tic, e.sched.tasks[l].toc,
                  (e.sched.tasks[l].ci != NULL) ? 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,
                  e.sched.tasks[l].flags);
            }
            fflush(stdout);
            count++;
          }
          fclose(file_thread);
        }

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

635
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
      char dumpfile[30];
      snprintf(dumpfile, 30, "thread_info-step%d.dat", j + 1);
      FILE *file_thread;
      file_thread = fopen(dumpfile, "w");
      /* Add some information to help with the plots */
      fprintf(file_thread, " %i %i %i %i %lli %lli %i %i %i %lli\n", -2, -1, -1,
              1, e.tic_step, e.toc_step, 0, 0, 0, cpufreq);
      for (int l = 0; l < e.sched.nr_tasks; l++) {
        if (!e.sched.tasks[l].implicit && e.sched.tasks[l].toc != 0) {
          fprintf(
              file_thread, " %i %i %i %i %lli %lli %i %i %i %i\n",
              e.sched.tasks[l].rid, e.sched.tasks[l].type,
              e.sched.tasks[l].subtype, (e.sched.tasks[l].cj == NULL),
              e.sched.tasks[l].tic, e.sched.tasks[l].toc,
              (e.sched.tasks[l].ci == NULL) ? 0 : e.sched.tasks[l].ci->count,
              (e.sched.tasks[l].cj == NULL) ? 0 : e.sched.tasks[l].cj->count,
              (e.sched.tasks[l].ci == NULL) ? 0 : e.sched.tasks[l].ci->gcount,
              (e.sched.tasks[l].cj == NULL) ? 0 : e.sched.tasks[l].cj->gcount);
        }
655
      }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
656
      fclose(file_thread);
657
#endif  // WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
658
    }
659
#endif  // SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
660
  }
661
662

/* Print the values of the runner histogram. */
663
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
664
665
666
667
668
669
670
  printf("main: runner histogram data:\n");
  for (k = 0; k < runner_hist_N; k++)
    printf(" %e %e %e\n",
           runner_hist_a + k * (runner_hist_b - runner_hist_a) / runner_hist_N,
           runner_hist_a +
               (k + 1) * (runner_hist_b - runner_hist_a) / runner_hist_N,
           (double)runner_hist_bins[k]);
671
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
672

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
673
  /* Write final output. */
674
  engine_drift_all(&e);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
675
  engine_dump_snapshot(&e);
676

677
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
678
679
  if ((res = MPI_Finalize()) != MPI_SUCCESS)
    error("call to MPI_Finalize failed with error %i.", res);
680
#endif
681

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
682
683
684
  /* Clean everything */
  engine_clean(&e);
  free(params);
685

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
689
690
  /* All is calm, all is bright. */
  return 0;
691
}