main.c 26 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.");
85
86
  printf("  %2s %8s %s\n", "-M", "",
         "Reconstruct the multipoles every time-step.");
Peter W. Draper's avatar
Peter W. Draper committed
87
  printf("  %2s %8s %s\n", "-n", "{int}",
88
89
         "Execute a fixed number of time steps. When unset use the time_end "
         "parameter to stop.");
90
91
  printf("  %2s %8s %s\n", "-P", "{sec:par:val}",
         "Set parameter value, can be used more than once.");
92
93
  printf("  %2s %8s %s\n", "-s", "", "Run with hydrodynamics.");
  printf("  %2s %8s %s\n", "-S", "", "Run with stars.");
94
95
96
  printf("  %2s %8s %s\n", "-t", "{int}",
         "The number of threads to use on each MPI rank. Defaults to 1 if not "
         "specified.");
97
  printf("  %2s %8s %s\n", "-T", "", "Print timers every time-step.");
98
  printf("  %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity.");
Matthieu Schaller's avatar
Matthieu Schaller committed
99
  printf("  %2s %8s %s\n", "", "", "1: MPI-rank 0 writes ");
100
  printf("  %2s %8s %s\n", "", "", "2: All MPI-ranks write");
101
  printf("  %2s %8s %s\n", "-y", "{int}",
102
103
         "Time-step frequency at which task graphs are dumped.");
  printf("  %2s %8s %s\n", "-h", "", "Print this help message and exit.");
104
105
106
  printf(
      "\nSee the file parameter_example.yml for an example of "
      "parameter file.\n");
107
108
}

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

115
  struct clocks_time tic, toc;
116

117
  int nr_nodes = 1, myrank = 0;
118

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

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

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

147
#endif
148

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    /* 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);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
421
422
423
  if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
  fflush(stdout);

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

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

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

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

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

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

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

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

  /* Initialize the engine with the space and policies. */
  if (myrank == 0) clocks_gettime(&tic);
  struct engine e;
559
  engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, N_total[0],
Peter W. Draper's avatar
Peter W. Draper committed
560
561
562
              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
563
564
565
566
567
568
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
569

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

575
576
577
578
579
580
581
582
583
584
585
586
#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
587
588
  /* Get some info to the user. */
  if (myrank == 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
589
    long long N_DM = N_total[1] - N_total[2] - N_total[0];
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
590
    message(
Matthieu Schaller's avatar
Matthieu Schaller committed
591
592
593
594
595
596
597
598
        "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
599
600
    fflush(stdout);
  }
601

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

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

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

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

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

633
  /* File for the timers */
634
  if (with_verbose_timers) timers_open_file(myrank);
635

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

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

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

645
    /* Print the timers. */
646
    if (with_verbose_timers) timers_print(e.step);
647

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
653
654
655
656
657
658
      /* 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");
659
        fclose(file_thread);
660
      }
661
      MPI_Barrier(MPI_COMM_WORLD);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
662
663
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

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

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

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

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

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

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

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

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