main.c 24.3 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
  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
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
338 339 340 341 342 343 344 345 346 347 348 349 350
  struct partition initial_partition;
  enum repartition_type reparttype;
  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
  /* Initialize unit system and constants */
354
  struct unit_system us;
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
355 356 357 358 359 360 361 362 363 364 365 366 367 368
  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
  /* Read particles and space information from (GADGET) ICs */
  char ICfileName[200] = "";
  parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
374 375
  const int replicate =
      parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
376 377 378
  if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
  fflush(stdout);

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

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

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
425
  /* Get the total number of particles across all nodes. */
426
  long long N_total[3] = {0, 0, 0};
427
#if defined(WITH_MPI)
428
  long long N_long[3] = {Ngas, Ngpart, Nspart};
Matthieu Schaller's avatar
Matthieu Schaller committed
429 430
  MPI_Reduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM, 0,
             MPI_COMM_WORLD);
431
#else
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
432 433
  N_total[0] = Ngas;
  N_total[1] = Ngpart;
434
  N_total[2] = Nspart;
435
#endif
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
436
  if (myrank == 0)
437 438 439 440
    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
441 442 443 444

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

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
454 455 456 457 458 459 460 461 462
  /* 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);
463
    message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
464 465 466
    message("maximum depth is %d.", s.maxdepth);
    fflush(stdout);
  }
467

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

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

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
482 483 484
  /* Initialise the external potential properties */
  struct external_potential potential;
  if (with_external_gravity)
485
    potential_init(params, &prog_const, &us, &s, &potential);
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506
  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;
507
  if (with_stars) engine_policies |= engine_policy_stars;
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
508 509 510 511 512 513 514 515 516 517 518 519 520

  /* 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,
              engine_policies, talking, &us, &prog_const, &hydro_properties,
              &potential, &cooling_func, &sourceterms);
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
521

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

527 528 529 530 531 532 533 534 535 536 537 538
#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
539 540
  /* Get some info to the user. */
  if (myrank == 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
541
    long long N_DM = N_total[1] - N_total[2] - N_total[0];
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
542
    message(
Matthieu Schaller's avatar
Matthieu Schaller committed
543 544 545 546 547 548 549 550
        "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
551 552
    fflush(stdout);
  }
553

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
554 555
  /* Time to say good-bye if this was not a serious run. */
  if (dry_run) {
556
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
557 558
    if ((res = MPI_Finalize()) != MPI_SUCCESS)
      error("call to MPI_Finalize failed with error %i.", res);
559
#endif
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
560 561 562 563 564 565
    if (myrank == 0)
      message("Time integration ready to start. End of dry-run.");
    engine_clean(&e);
    free(params);
    return 0;
  }
566 567

#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
568 569 570
  /* Split the space. */
  engine_split(&e, &initial_partition);
  engine_redistribute(&e);
571 572
#endif

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
573 574
  /* Initialise the particles */
  engine_init_particles(&e, flag_entropy_ICs);
575

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

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
579 580
  /* Legend */
  if (myrank == 0)
581 582 583
    printf("# %6s %14s %14s %10s %10s %10s %16s [%s]\n", "Step", "Time",
           "Time-step", "Updates", "g-Updates", "s-Updates", "Wall-clock time",
           clocks_getunit());
584

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

/* Repartition the space amongst the nodes? */
589
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
590
    if (j % 100 == 2) e.forcerepart = reparttype;
591 592
#endif

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
593 594
    /* Reset timers */
    timers_reset(timers_mask_all);
595

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
596 597
    /* Take a step. */
    engine_step(&e);
598

599
#ifdef SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
600 601
    /* Dump the task data using the given frequency. */
    if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
602 603
#ifdef WITH_MPI

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
604 605 606 607 608 609
      /* 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");
610
        fclose(file_thread);
611
      }
612
      MPI_Barrier(MPI_COMM_WORLD);
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653

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

655
#else
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674
      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);
        }
675
      }
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
676
      fclose(file_thread);
677
#endif  // WITH_MPI
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
678
    }
679
#endif  // SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
680
  }
681 682

/* Print the values of the runner histogram. */
683
#ifdef HIST
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
684 685 686 687 688 689 690
  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]);
691
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
692

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
693 694
  /* Write final output. */
  engine_dump_snapshot(&e);
695

696
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
697 698
  if ((res = MPI_Finalize()) != MPI_SUCCESS)
    error("call to MPI_Finalize failed with error %i.", res);
699
#endif
700

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
701 702 703
  /* Clean everything */
  engine_clean(&e);
  free(params);
704

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
705 706
  /* Say goodbye. */
  if (myrank == 0) message("done. Bye.");
707

Pedro Gonnet's avatar
oops.  
Pedro Gonnet committed
708 709
  /* All is calm, all is bright. */
  return 0;
710
}