main.c 26.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 <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
  printf("\nUsage: swift [OPTION]... PARAMFILE\n");
57
  printf("       swift_mpi [OPTION]... PARAMFILE\n\n");
58

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

Pedro Gonnet's avatar
Pedro Gonnet committed
110
111
112
113
/**
 * @brief Main routine that loads a few particles and generates some output.
 *
 */
114
115
int main(int argc, char *argv[]) {

116
  struct clocks_time tic, toc;
117

118
  int nr_nodes = 1, myrank = 0;
119

120
#ifdef WITH_MPI
121
  /* Start by initializing MPI. */
122
  int res = 0, prov = 0;
123
124
  if ((res = MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &prov)) !=
      MPI_SUCCESS)
125
126
    error("Call to MPI_Init failed with error %i.", res);
  if (prov != MPI_THREAD_MULTIPLE)
127
    error(
128
129
        "MPI does not provide the level of threading"
        " required (MPI_THREAD_MULTIPLE).");
130
  if ((res = MPI_Comm_size(MPI_COMM_WORLD, &nr_nodes)) != MPI_SUCCESS)
131
132
133
    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);
134
135
136
137

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

138
139
140
  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);
141
  if (myrank == 0) message("MPI is up and running with %i node(s).", nr_nodes);
142
143
144
145
  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.");
  }
146
  fflush(stdout);
147

148
#endif
149

150
/* Let's pin the main thread */
151
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
152
  if (((ENGINE_POLICY)&engine_policy_setaffinity) == engine_policy_setaffinity)
153
    engine_pin();
Angus Lepper's avatar
Angus Lepper committed
154
155
#endif

156
157
  /* Welcome to SWIFT, you made the right choice */
  if (myrank == 0) greetings();
158

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

  /* Parse the parameters */
182
  int c;
183
  while ((c = getopt(argc, argv, "acCdDef:FgGhMn:P:sSt:Tv:y:")) != -1)
184
    switch (c) {
185
      case 'a':
186
        with_aff = 1;
187
        break;
188
      case 'c':
189
        with_cosmology = 1;
190
        break;
191
192
193
      case 'C':
        with_cooling = 1;
        break;
194
      case 'd':
195
        dry_run = 1;
196
        break;
197
198
199
      case 'D':
        with_drift_all = 1;
        break;
200
      case 'e':
201
        with_fp_exceptions = 1;
202
        break;
203
      case 'f':
204
205
206
        if (sscanf(optarg, "%llu", &cpufreq) != 1) {
          if (myrank == 0) printf("Error parsing CPU frequency (-f).\n");
          if (myrank == 0) print_help_message();
207
          return 1;
208
        }
209
        break;
210
211
212
      case 'F':
        with_sourceterms = 1;
        break;
213
      case 'g':
214
        with_external_gravity = 1;
215
        break;
216
217
      case 'G':
        with_self_gravity = 1;
Tom Theuns's avatar
Tom Theuns committed
218
        break;
219
      case 'h':
220
        if (myrank == 0) print_help_message();
221
        return 0;
222
223
224
      case 'M':
        with_mpole_reconstruction = 1;
        break;
225
226
227
228
229
230
231
      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;
232
233
234
235
      case 'P':
        cmdparams[nparams] = optarg;
        nparams++;
        break;
236
      case 's':
237
        with_hydro = 1;
238
        break;
239
240
241
      case 'S':
        with_stars = 1;
        break;
242
243
244
245
246
247
248
249
      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;
250
251
252
      case 'T':
        with_verbose_timers = 1;
        break;
253
      case 'v':
254
255
256
        if (sscanf(optarg, "%d", &verbose) != 1) {
          if (myrank == 0) printf("Error parsing verbosity level (-v).\n");
          if (myrank == 0) print_help_message();
257
          return 1;
258
        }
259
        break;
260
      case 'y':
261
262
263
        if (sscanf(optarg, "%d", &dump_tasks) != 1) {
          if (myrank == 0) printf("Error parsing dump_tasks (-y). \n");
          if (myrank == 0) print_help_message();
264
          return 1;
265
        }
266
267
268
269
270
271
272
#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
273
274
275
276
277
        break;
      case '?':
        if (myrank == 0) print_help_message();
        return 1;
        break;
278
    }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
279
280
281
282
283
  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");
284
285
    if (myrank == 0) print_help_message();
    return 1;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
286
287
288
289
290
291
292
293
294
295
296
  } 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;
  }
297

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

301
302
303
  /* How vocal are we ? */
  const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
308
309
310
311
312
  /* Report CPU frequency.*/
  cpufreq = clocks_get_cpufreq();
  if (myrank == 0) {
    message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
  }
313

Matthieu Schaller's avatar
Matthieu Schaller committed
314
/* Report host name(s). */
315
#ifdef WITH_MPI
316
  if (talking) {
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
317
318
    message("Rank %d running on: %s", myrank, hostname());
  }
319
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
320
  message("Running on: %s", hostname());
321
322
#endif

323
324
/* Do we have debugging checks ? */
#ifdef SWIFT_DEBUG_CHECKS
325
326
  if (myrank == 0)
    message("WARNING: Debugging checks activated. Code will be slower !");
327
328
#endif

329
330
331
332
333
334
335
336
337
/* 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
338
339
340
  /* Do we choke on FP-exceptions ? */
  if (with_fp_exceptions) {
    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
341
342
    if (myrank == 0)
      message("WARNING: Floating point exceptions will be reported.");
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
343
  }
344

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
345
346
  /* How large are the parts? */
  if (myrank == 0) {
347
348
349
350
351
352
353
354
    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
355
  }
356

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
357
358
359
360
361
362
  /* 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);
363
364
365

    /* Handle any command-line overrides. */
    if (nparams > 0)
366
      for (int k = 0; k < nparams; k++) parser_set_param(params, cmdparams[k]);
367
368

    /* And dump the parameters as used. */
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
369
370
371
    // parser_print_params(&params);
    parser_write_params_to_file(params, "used_parameters.yml");
  }
372
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
373
374
  /* Broadcast the parameter file */
  MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
375
#endif
376

377
  /* Prepare the domain decomposition scheme */
378
  struct repartition reparttype;
379
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
380
381
382
383
384
385
386
387
388
389
  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]);
390
    message("Using %s repartitioning", repartition_name[reparttype.type]);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
391
  }
392
#endif
393

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
394
  /* Initialize unit system and constants */
395
  struct unit_system us;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
396
397
398
399
400
401
402
403
404
405
406
407
408
409
  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
410
  if (with_hydro) hydro_props_init(&hydro_properties, params);
411

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

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

426
  /* Get ready to read particles of all kinds */
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
427
428
  struct part *parts = NULL;
  struct gpart *gparts = NULL;
429
430
  struct spart *sparts = NULL;
  size_t Ngas = 0, Ngpart = 0, Nspart = 0;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
431
432
433
434
  double dim[3] = {0., 0., 0.};
  int periodic = 0;
  int flag_entropy_ICs = 0;
  if (myrank == 0) clocks_gettime(&tic);
435
436
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
Matthieu Schaller's avatar
Matthieu Schaller committed
437
438
439
440
  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
441
#else
442
  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
Matthieu Schaller's avatar
Typos    
Matthieu Schaller committed
443
                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
444
445
                 (with_external_gravity || with_self_gravity), with_stars,
                 myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
446
447
#endif
#else
448
  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
449
450
451
                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                 (with_external_gravity || with_self_gravity), with_stars,
                 dry_run);
452
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
453
454
455
456
457
458
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("Reading initial conditions took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
459

460
461
#ifdef SWIFT_DEBUG_CHECKS
  /* Check once and for all that we don't have unwanted links */
462
  if (!with_stars) {
463
    for (size_t k = 0; k < Ngpart; ++k)
464
465
      if (gparts[k].type == swift_type_star) error("Linking problem");
  }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
466
467
  if (!with_hydro) {
    for (size_t k = 0; k < Ngpart; ++k)
468
      if (gparts[k].type == swift_type_gas) error("Linking problem");
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
469
  }
470
#endif
471

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
472
  /* Get the total number of particles across all nodes. */
473
  long long N_total[3] = {0, 0, 0};
474
#if defined(WITH_MPI)
475
  long long N_long[3] = {Ngas, Ngpart, Nspart};
476
477
  MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
                MPI_COMM_WORLD);
478
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
479
480
  N_total[0] = Ngas;
  N_total[1] = Ngpart;
481
  N_total[2] = Nspart;
482
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
483
  if (myrank == 0)
484
485
486
487
    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
488
489
490
491

  /* Initialize the space with these data. */
  if (myrank == 0) clocks_gettime(&tic);
  struct space s;
492
  space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
493
             periodic, replicate, with_self_gravity, talking, dry_run);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
494
495
496
497
498
499
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
500

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
501
502
503
504
505
506
507
508
509
  /* 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);
510
    message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
511
512
513
    message("maximum depth is %d.", s.maxdepth);
    fflush(stdout);
  }
514

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
515
516
517
518
519
520
  /* 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);
  }
521

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
522
523
524
525
526
527
  /* 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]);
  }
528

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
529
530
531
  /* Initialise the external potential properties */
  struct external_potential potential;
  if (with_external_gravity)
532
    potential_init(params, &prog_const, &us, &s, &potential);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
  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;
548
549
  if (with_mpole_reconstruction)
    engine_policies |= engine_policy_reconstruct_mpoles;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
550
551
552
553
554
555
  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;
556
  if (with_stars) engine_policies |= engine_policy_stars;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
557
558
559
560

  /* Initialize the engine with the space and policies. */
  if (myrank == 0) clocks_gettime(&tic);
  struct engine e;
561
  engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, N_total[0],
Peter W. Draper's avatar
Peter W. Draper committed
562
563
564
              N_total[1], with_aff, engine_policies, talking, &reparttype, &us,
              &prog_const, &hydro_properties, &gravity_properties, &potential,
              &cooling_func, &sourceterms);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
565
566
567
568
569
570
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
571

572
573
/* Init the runner history. */
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
574
  for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
575
576
#endif

577
578
579
580
581
582
583
584
585
586
587
588
#if defined(WITH_MPI)
  N_long[0] = s.nr_parts;
  N_long[1] = s.nr_gparts;
  N_long[2] = s.nr_sparts;
  MPI_Reduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM, 0,
             MPI_COMM_WORLD);
#else
  N_total[0] = s.nr_parts;
  N_total[1] = s.nr_gparts;
  N_total[2] = s.nr_sparts;
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
589
590
  /* Get some info to the user. */
  if (myrank == 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
591
    long long N_DM = N_total[1] - N_total[2] - N_total[0];
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
592
    message(
Matthieu Schaller's avatar
Matthieu Schaller committed
593
594
595
596
597
598
599
600
        "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
601
602
    fflush(stdout);
  }
603

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
604
605
  /* Time to say good-bye if this was not a serious run. */
  if (dry_run) {
606
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
607
608
    if ((res = MPI_Finalize()) != MPI_SUCCESS)
      error("call to MPI_Finalize failed with error %i.", res);
609
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
610
611
612
613
614
615
    if (myrank == 0)
      message("Time integration ready to start. End of dry-run.");
    engine_clean(&e);
    free(params);
    return 0;
  }
616
617

#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
618
619
620
  /* Split the space. */
  engine_split(&e, &initial_partition);
  engine_redistribute(&e);
621
622
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
623
  /* Initialise the particles */
624
  engine_init_particles(&e, flag_entropy_ICs, clean_h_values);
625

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
629
  /* Legend */
630
  if (myrank == 0)
631
632
633
    printf("# %6s %14s %14s %10s %10s %10s %16s [%s]\n", "Step", "Time",
           "Time-step", "Updates", "g-Updates", "s-Updates", "Wall-clock time",
           clocks_getunit());
634

635
  /* File for the timers */
636
  if (with_verbose_timers) timers_open_file(myrank);
637

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
641
    /* Reset timers */
642
    timers_reset_all();
643

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
644
645
    /* Take a step. */
    engine_step(&e);
646

647
    /* Print the timers. */
648
    if (with_verbose_timers) timers_print(e.step);
649

650
#ifdef SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
651
652
    /* Dump the task data using the given frequency. */
    if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
653
654
#ifdef WITH_MPI

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
655
656
657
658
659
660
      /* 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");
661
        fclose(file_thread);
662
      }
663
      MPI_Barrier(MPI_COMM_WORLD);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704

      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);
      }
705

706
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
      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);
        }
726
      }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
727
      fclose(file_thread);
728
#endif  // WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
729
    }
730
#endif  // SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
731
  }
732
733

/* Print the values of the runner histogram. */
734
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
735
736
737
738
739
740
741
  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]);
742
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
743

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
744
  /* Write final output. */
745
  engine_drift_all(&e);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
746
  engine_dump_snapshot(&e);
747

748
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
749
750
  if ((res = MPI_Finalize()) != MPI_SUCCESS)
    error("call to MPI_Finalize failed with error %i.", res);
751
#endif
752

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
753
  /* Clean everything */
754
  if (with_verbose_timers) timers_close_file();
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
755
756
  engine_clean(&e);
  free(params);
757

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
761
762
  /* All is calm, all is bright. */
  return 0;
763
}