main.c 25.8 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
  printf("  %2s %8s %s\n", "-a", "", "Pin runners using processor affinity.");
  printf("  %2s %8s %s\n", "-c", "", "Run with cosmological time integration.");
  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
  printf("  %2s %8s %s\n", "-D", "",
73
74
75
76
         "Always drift all particles even the ones far from active particles. "
         "This emulates");
  printf("  %2s %8s %s\n", "", "",
         "Gadget-[23] and GIZMO's default behaviours.");
77
  printf("  %2s %8s %s\n", "-e", "",
78
         "Enable floating-point exceptions (debugging mode).");
79
  printf("  %2s %8s %s\n", "-f", "{int}",
80
         "Overwrite the CPU frequency (Hz) to be used for time measurements.");
81
  printf("  %2s %8s %s\n", "-g", "",
82
83
84
         "Run with an external gravitational potential.");
  printf("  %2s %8s %s\n", "-F", "", "Run with feedback.");
  printf("  %2s %8s %s\n", "-G", "", "Run with self-gravity.");
Peter W. Draper's avatar
Peter W. Draper committed
85
  printf("  %2s %8s %s\n", "-n", "{int}",
86
87
         "Execute a fixed number of time steps. When unset use the time_end "
         "parameter to stop.");
88
89
  printf("  %2s %8s %s\n", "-s", "", "Run with hydrodynamics.");
  printf("  %2s %8s %s\n", "-S", "", "Run with stars.");
90
91
92
  printf("  %2s %8s %s\n", "-t", "{int}",
         "The number of threads to use on each MPI rank. Defaults to 1 if not "
         "specified.");
93
  printf("  %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity.");
Matthieu Schaller's avatar
Matthieu Schaller committed
94
  printf("  %2s %8s %s\n", "", "", "1: MPI-rank 0 writes ");
95
  printf("  %2s %8s %s\n", "", "", "2: All MPI-ranks write");
96
  printf("  %2s %8s %s\n", "-y", "{int}",
97
98
         "Time-step frequency at which task graphs are dumped.");
  printf("  %2s %8s %s\n", "-h", "", "Print this help message and exit.");
99
100
101
  printf(
      "\nSee the file parameter_example.yml for an example of "
      "parameter file.\n");
102
103
}

104
105
106
107
#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
VORONOI3D_DECLARE_GLOBAL_VARIABLES()
#endif

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

114
  struct clocks_time tic, toc;
115

116
  int nr_nodes = 1, myrank = 0;
117

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

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

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

146
#endif
147

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

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

157
  int with_aff = 0;
158
  int dry_run = 0;
159
  int dump_tasks = 0;
160
  int nsteps = -2;
161
162
  int with_cosmology = 0;
  int with_external_gravity = 0;
Tom Theuns's avatar
Tom Theuns committed
163
  int with_sourceterms = 0;
164
  int with_cooling = 0;
165
166
  int with_self_gravity = 0;
  int with_hydro = 0;
167
  int with_stars = 0;
168
  int with_fp_exceptions = 0;
169
  int with_drift_all = 0;
170
  int verbose = 0;
171
  int nr_threads = 1;
172
173
174
175
  char paramFileName[200] = "";
  unsigned long long cpufreq = 0;

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

281
  /* Genesis 1.1: And then, there was time ! */
282
  clocks_set_cpufreq(cpufreq);
283

284
285
286
  /* How vocal are we ? */
  const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);

287
288
289
290
  if (myrank == 0 && dry_run)
    message(
        "Executing a dry run. No i/o or time integration will be performed.");

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
291
  /* Report CPU frequency.*/
292
  cpufreq = clocks_get_cpufreq();
293
  if (myrank == 0) {
294
    message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
295
296
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
297
/* Report host name(s). */
298
#ifdef WITH_MPI
299
  if (talking) {
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
300
301
    message("Rank %d running on: %s", myrank, hostname());
  }
302
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
303
  message("Running on: %s", hostname());
304
305
#endif

306
307
/* Do we have debugging checks ? */
#ifdef SWIFT_DEBUG_CHECKS
308
309
  if (myrank == 0)
    message("WARNING: Debugging checks activated. Code will be slower !");
310
311
#endif

312
313
314
315
316
317
318
319
320
/* 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

321
322
  /* Do we choke on FP-exceptions ? */
  if (with_fp_exceptions) {
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
323
    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
324
325
    if (myrank == 0)
      message("WARNING: Floating point exceptions will be reported.");
326
  }
327

328
329
  /* How large are the parts? */
  if (myrank == 0) {
330
331
332
333
334
335
336
337
    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(acc_tensor) is %4zi bytes.", sizeof(struct acc_tensor));
    message("sizeof(task)       is %4zi bytes.", sizeof(struct task));
    message("sizeof(cell)       is %4zi bytes.", sizeof(struct cell));
338
  }
339

340
  /* Read the parameter file */
341
342
  struct swift_params *params = malloc(sizeof(struct swift_params));
  if (params == NULL) error("Error allocating memory for the parameter file.");
343
  if (myrank == 0) {
344
    message("Reading runtime parameters from file '%s'", paramFileName);
345
    parser_read_file(paramFileName, params);
346
    // parser_print_params(&params);
347
    parser_write_params_to_file(params, "used_parameters.yml");
348
  }
349
350
#ifdef WITH_MPI
  /* Broadcast the parameter file */
351
  MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
352
#endif
353

354
  /* Prepare the domain decomposition scheme */
355
  enum repartition_type reparttype = REPART_NONE;
356
#ifdef WITH_MPI
357
  struct partition initial_partition;
358
  partition_init(&initial_partition, &reparttype, params, nr_nodes);
359
360

  /* Let's report what we did */
361
  if (myrank == 0) {
362
363
364
365
366
367
    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]);
368
  }
369
#endif
370

Matthieu Schaller's avatar
Matthieu Schaller committed
371
  /* Initialize unit system and constants */
372
  struct unit_system us;
373
  struct phys_const prog_const;
374
  units_init(&us, params, "InternalUnitSystem");
375
376
  phys_const_init(&us, &prog_const);
  if (myrank == 0 && verbose > 0) {
377
378
379
380
381
    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);
382
    phys_const_print(&prog_const);
383
  }
384

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

389
390
391
  /* Initialise the gravity properties */
  struct gravity_props gravity_properties;
  if (with_self_gravity) gravity_props_init(&gravity_properties, params);
392

393
#if defined(SHADOWFAX_SPH)
394
395
396
397
398
399
  /* Override the variables governing the density iteration
     (see src/hydro/Shadowswift/hydro.h for full explanation) */
  hydro_properties.target_neighbours = 1.0f;
  hydro_properties.delta_neighbours = 0.0f;
#endif

400
  /* Read particles and space information from (GADGET) ICs */
401
  char ICfileName[200] = "";
402
  parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
403
404
  const int replicate =
      parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
405
  if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
406
407
  fflush(stdout);

408
  /* Get ready to read particles of all kinds */
409
410
  struct part *parts = NULL;
  struct gpart *gparts = NULL;
411
412
  struct spart *sparts = NULL;
  size_t Ngas = 0, Ngpart = 0, Nspart = 0;
413
414
  double dim[3] = {0., 0., 0.};
  int periodic = 0;
415
  int flag_entropy_ICs = 0;
416
  if (myrank == 0) clocks_gettime(&tic);
417
418
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
Matthieu Schaller's avatar
Matthieu Schaller committed
419
420
421
422
  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);
423
#else
424
  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
Matthieu Schaller's avatar
Typos    
Matthieu Schaller committed
425
                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
426
427
                 (with_external_gravity || with_self_gravity), with_stars,
                 myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
428
429
#endif
#else
430
  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
431
432
433
                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                 (with_external_gravity || with_self_gravity), with_stars,
                 dry_run);
434
#endif
435
436
  if (myrank == 0) {
    clocks_gettime(&toc);
437
438
    message("Reading initial conditions took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
439
440
    fflush(stdout);
  }
441

442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
  /* set the *global* box dimensions */
  float box_anchor[3], box_side[3];
  if (periodic) {
    box_anchor[0] = -0.5f * dim[0];
    box_anchor[1] = -0.5f * dim[1];
    box_anchor[2] = -0.5f * dim[2];
    box_side[0] = 2.0f * dim[0];
    box_side[1] = 2.0f * dim[1];
    box_side[2] = 2.0f * dim[2];
  } else {
    box_anchor[0] = 0.0f;
    box_anchor[1] = 0.0f;
    box_anchor[2] = 0.0f;
    box_side[0] = dim[0];
    box_side[1] = dim[1];
    box_side[2] = dim[2];
  }
  voronoi_set_box(box_anchor, box_side);
#endif

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

  /* Get the total number of particles across all nodes. */
476
  long long N_total[3] = {0, 0, 0};
477
#if defined(WITH_MPI)
478
  long long N_long[3] = {Ngas, Ngpart, Nspart};
Matthieu Schaller's avatar
Matthieu Schaller committed
479
480
  MPI_Reduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM, 0,
             MPI_COMM_WORLD);
481
#else
482
  N_total[0] = Ngas;
483
  N_total[1] = Ngpart;
484
  N_total[2] = Nspart;
485
#endif
486
  if (myrank == 0)
487
488
489
490
    message(
        "Read %lld gas particles, %lld star particles and %lld gparts from the "
        "ICs.",
        N_total[0], N_total[2], N_total[1]);
491

492
  /* Initialize the space with these data. */
493
  if (myrank == 0) clocks_gettime(&tic);
494
  struct space s;
495
  space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
496
             periodic, replicate, with_self_gravity, talking, dry_run);
497
  if (myrank == 0) {
498
    clocks_gettime(&toc);
499
500
    message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
501
502
    fflush(stdout);
  }
503
504
505
506
507
508
509
510

  /* 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]);
Peter W. Draper's avatar
Peter W. Draper committed
511
    message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
512
    message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
513
    message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
Matthieu Schaller's avatar
Matthieu Schaller committed
514
515
    message("maximum depth is %d.", s.maxdepth);
    fflush(stdout);
516
517
  }

518
  /* Verify that each particle is in it's proper cell. */
519
  if (talking && !dry_run) {
520
    int icount = 0;
521
522
523
524
    space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
    message("map_cellcheck picked up %i parts.", icount);
  }

525
  /* Verify the maximal depth of cells. */
526
  if (talking && !dry_run) {
527
    int data[2] = {s.maxdepth, 0};
528
529
530
531
    space_map_cells_pre(&s, 0, &map_maxdepth, data);
    message("nr of cells at depth %i is %i.", data[0], data[1]);
  }

532
533
534
  /* Initialise the external potential properties */
  struct external_potential potential;
  if (with_external_gravity)
535
    potential_init(params, &prog_const, &us, &s, &potential);
536
537
538
  if (with_external_gravity && myrank == 0) potential_print(&potential);

  /* Initialise the cooling function properties */
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
539
540
541
542
543
544
545
546
  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);
547

548
  /* Construct the engine policy */
549
  int engine_policies = ENGINE_POLICY | engine_policy_steal;
550
  if (with_drift_all) engine_policies |= engine_policy_drift_all;
551
  if (with_hydro) engine_policies |= engine_policy_hydro;
552
  if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
553
554
  if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
  if (with_cosmology) engine_policies |= engine_policy_cosmology;
555
  if (with_cooling) engine_policies |= engine_policy_cooling;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
556
  if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
557
  if (with_stars) engine_policies |= engine_policy_stars;
558

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

573
574
575
576
577
/* Init the runner history. */
#ifdef HIST
  for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
#endif

578
579
580
581
582
583
584
585
586
587
588
589
#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

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

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

#ifdef WITH_MPI
  /* Split the space. */
  engine_split(&e, &initial_partition);
  engine_redistribute(&e);
#endif

624
  /* Initialise the particles */
625
  engine_init_particles(&e, flag_entropy_ICs);
626

627
628
629
  /* Write the state of the system before starting time integration. */
  engine_dump_snapshot(&e);

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

636
  /* Main simulation loop */
637
  for (int j = 0; !engine_is_done(&e) && e.step - 1 != nsteps; j++) {
638

639
    /* Reset timers */
640
641
642
643
644
    timers_reset(timers_mask_all);

    /* Take a step. */
    engine_step(&e);

645
#ifdef SWIFT_DEBUG_TASKS
646
    /* Dump the task data using the given frequency. */
647
648
649
    if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
#ifdef WITH_MPI

650
      /* Make sure output file is empty, only on one rank. */
651
      char dumpfile[30];
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
652
      snprintf(dumpfile, 30, "thread_info_MPI-step%d.dat", j + 1);
653
      FILE *file_thread;
654
      if (myrank == 0) {
655
656
        file_thread = fopen(dumpfile, "w");
        fclose(file_thread);
657
658
      }
      MPI_Barrier(MPI_COMM_WORLD);
659
660
661
662
663
664
665
666
667
668
669
670

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

671
672
          fprintf(file_thread, " %03i 0 0 0 0 %lli %lli 0 0 0 0 %lli\n", myrank,
                  e.tic_step, e.toc_step, cpufreq);
673
          int count = 0;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
674
675
          for (int l = 0; l < e.sched.nr_tasks; l++) {
            if (!e.sched.tasks[l].implicit && e.sched.tasks[l].toc != 0) {
676
677
              fprintf(
                  file_thread, " %03i %i %i %i %i %lli %lli %i %i %i %i %i\n",
678
                  myrank, e.sched.tasks[l].rid, e.sched.tasks[l].type,
679
680
681
682
683
684
685
686
687
688
689
                  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);
690
            }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
691
692
693
            fflush(stdout);
            count++;
          }
694
695
696
697
698
699
700
701
          fclose(file_thread);
        }

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

#else
702
      char dumpfile[30];
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
703
      snprintf(dumpfile, 30, "thread_info-step%d.dat", j + 1);
704
      FILE *file_thread;
705
      file_thread = fopen(dumpfile, "w");
706
      /* Add some information to help with the plots */
707
708
      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);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
709
710
      for (int l = 0; l < e.sched.nr_tasks; l++) {
        if (!e.sched.tasks[l].implicit && e.sched.tasks[l].toc != 0) {
711
          fprintf(
712
              file_thread, " %i %i %i %i %lli %lli %i %i %i %i\n",
713
              e.sched.tasks[l].rid, e.sched.tasks[l].type,
714
715
716
              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,
717
718
719
              (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);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
720
        }
721
      }
722
      fclose(file_thread);
723
#endif  // WITH_MPI
724
    }
725
#endif  // SWIFT_DEBUG_TASKS
726
727
728
  }

/* Print the values of the runner histogram. */
729
#ifdef HIST
730
731
732
733
734
735
736
  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]);
737
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
738

739
  /* Write final output. */
740
  engine_drift_all(&e);
741
  engine_dump_snapshot(&e);
742

743
#ifdef WITH_MPI
744
  if ((res = MPI_Finalize()) != MPI_SUCCESS)
745
    error("call to MPI_Finalize failed with error %i.", res);
746
#endif
747

748
749
750
751
  /* Clean everything */
  engine_clean(&e);
  free(params);

752
  /* Say goodbye. */
Matthieu Schaller's avatar
Matthieu Schaller committed
753
  if (myrank == 0) message("done. Bye.");
754
755
756

  /* All is calm, all is bright. */
  return 0;
757
}