main.c 25.2 KB
Newer Older
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
4
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
5
 *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
6
7
8
 *                    Angus Lepper (angus.lepper@ed.ac.uk)
 *               2016 John A. Regan (john.a.regan@durham.ac.uk)
 *                    Tom Theuns (tom.theuns@durham.ac.uk)
9
 *
10
11
12
13
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
14
 *
15
16
17
18
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
19
 *
20
21
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
22
 *
23
24
 ******************************************************************************/

Pedro Gonnet's avatar
Pedro Gonnet committed
25
26
/* Config parameters. */
#include "../config.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
27
28

/* Some standard headers. */
29
#include <fenv.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
30
31
32
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
33
#include <unistd.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
34

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

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

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

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

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

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

59
  printf("Valid options are:\n");
60
61
62
  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
}

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

110
  struct clocks_time tic, toc;
111

112
  int nr_nodes = 1, myrank = 0;
113

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

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

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

142
#endif
143

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

150
151
  /* Welcome to SWIFT, you made the right choice */
  if (myrank == 0) greetings();
152

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

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

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

280
281
282
  /* How vocal are we ? */
  const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);

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

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

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

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

308
309
310
311
312
313
314
315
316
/* 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
317
318
319
  /* Do we choke on FP-exceptions ? */
  if (with_fp_exceptions) {
    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
320
321
    if (myrank == 0)
      message("WARNING: Floating point exceptions will be reported.");
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
322
  }
323

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
324
325
  /* How large are the parts? */
  if (myrank == 0) {
326
327
328
329
330
331
332
333
    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));
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
334
  }
335

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
336
337
338
339
340
341
342
343
344
  /* 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");
  }
345
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
346
347
  /* Broadcast the parameter file */
  MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
348
#endif
349

350
  /* Prepare the domain decomposition scheme */
351
  enum repartition_type reparttype = REPART_NONE;
352
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
353
354
355
356
357
358
359
360
361
362
363
364
  struct partition initial_partition;
  partition_init(&initial_partition, &reparttype, params, nr_nodes);

  /* Let's report what we did */
  if (myrank == 0) {
    message("Using initial partition %s",
            initial_partition_name[initial_partition.type]);
    if (initial_partition.type == INITPART_GRID)
      message("grid set to [ %i %i %i ].", initial_partition.grid[0],
              initial_partition.grid[1], initial_partition.grid[2]);
    message("Using %s repartitioning", repartition_name[reparttype]);
  }
365
#endif
366

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
367
  /* Initialize unit system and constants */
368
  struct unit_system us;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
369
370
371
372
373
374
375
376
377
378
379
380
381
382
  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
383
  if (with_hydro) hydro_props_init(&hydro_properties, params);
384

385
386
387
388
  /* 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
389
390
391
  /* Read particles and space information from (GADGET) ICs */
  char ICfileName[200] = "";
  parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
392
393
  const int replicate =
      parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
394
395
396
  if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
  fflush(stdout);

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

431
432
#ifdef SWIFT_DEBUG_CHECKS
  /* Check once and for all that we don't have unwanted links */
433
  if (!with_stars) {
434
    for (size_t k = 0; k < Ngpart; ++k)
435
436
      if (gparts[k].type == swift_type_star) error("Linking problem");
  }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
437
438
  if (!with_hydro) {
    for (size_t k = 0; k < Ngpart; ++k)
439
      if (gparts[k].type == swift_type_gas) error("Linking problem");
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
440
  }
441
#endif
442

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
443
  /* Get the total number of particles across all nodes. */
444
  long long N_total[3] = {0, 0, 0};
445
#if defined(WITH_MPI)
446
  long long N_long[3] = {Ngas, Ngpart, Nspart};
Matthieu Schaller's avatar
Matthieu Schaller committed
447
448
  MPI_Reduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM, 0,
             MPI_COMM_WORLD);
449
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
450
451
  N_total[0] = Ngas;
  N_total[1] = Ngpart;
452
  N_total[2] = Nspart;
453
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
454
  if (myrank == 0)
455
456
457
458
    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
459
460
461
462

  /* Initialize the space with these data. */
  if (myrank == 0) clocks_gettime(&tic);
  struct space s;
463
  space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
464
             periodic, replicate, with_self_gravity, talking, dry_run);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
465
466
467
468
469
470
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
471

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
472
473
474
475
476
477
478
479
480
  /* 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);
481
    message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
482
483
484
    message("maximum depth is %d.", s.maxdepth);
    fflush(stdout);
  }
485

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
486
487
488
489
490
491
  /* 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);
  }
492

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
493
494
495
496
497
498
  /* 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]);
  }
499

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
500
501
502
  /* Initialise the external potential properties */
  struct external_potential potential;
  if (with_external_gravity)
503
    potential_init(params, &prog_const, &us, &s, &potential);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
  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;
525
  if (with_stars) engine_policies |= engine_policy_stars;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
526
527
528
529
530

  /* 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,
531
              engine_policies, talking, reparttype, &us, &prog_const,
Matthieu Schaller's avatar
Matthieu Schaller committed
532
533
              &hydro_properties, &gravity_properties, &potential, &cooling_func,
              &sourceterms);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
534
535
536
537
538
539
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
540

541
542
/* Init the runner history. */
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
543
  for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
544
545
#endif

546
547
548
549
550
551
552
553
554
555
556
557
#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
558
559
  /* Get some info to the user. */
  if (myrank == 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
560
    long long N_DM = N_total[1] - N_total[2] - N_total[0];
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
561
    message(
Matthieu Schaller's avatar
Matthieu Schaller committed
562
563
564
565
566
567
568
569
        "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
570
571
    fflush(stdout);
  }
572

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
573
574
  /* Time to say good-bye if this was not a serious run. */
  if (dry_run) {
575
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
576
577
    if ((res = MPI_Finalize()) != MPI_SUCCESS)
      error("call to MPI_Finalize failed with error %i.", res);
578
#endif
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
579
580
581
582
583
584
    if (myrank == 0)
      message("Time integration ready to start. End of dry-run.");
    engine_clean(&e);
    free(params);
    return 0;
  }
585
586

#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
587
588
589
  /* Split the space. */
  engine_split(&e, &initial_partition);
  engine_redistribute(&e);
590
591
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
592
593
  /* Initialise the particles */
  engine_init_particles(&e, flag_entropy_ICs);
594

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
598
  /* Legend */
Pedro Gonnet's avatar
Pedro Gonnet committed
599
  if (myrank == 0) {
600
601
602
    printf("# %6s %14s %14s %10s %10s %10s %16s [%s]\n", "Step", "Time",
           "Time-step", "Updates", "g-Updates", "s-Updates", "Wall-clock time",
           clocks_getunit());
Pedro Gonnet's avatar
Pedro Gonnet committed
603
604
605
606
607
    printf("timers: ");
    for (int k = 0; k < timer_count; k++)
      printf("%s\t", timers_names[k]);
    printf("\n");
  }
608

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
612
613
    /* Reset timers */
    timers_reset(timers_mask_all);
614

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
615
616
    /* Take a step. */
    engine_step(&e);
Pedro Gonnet's avatar
Pedro Gonnet committed
617
618
619
620
621
622
    
    /* Print the timers. */
    printf("timers: ");
    for (int k = 0; k < timer_count; k++)
      printf("%.3f\t", clocks_from_ticks(timers[k]));
    printf("\n");
623

624
#ifdef SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
625
626
    /* Dump the task data using the given frequency. */
    if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
627
628
#ifdef WITH_MPI

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
629
630
631
632
633
634
      /* 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");
635
        fclose(file_thread);
636
      }
637
      MPI_Barrier(MPI_COMM_WORLD);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678

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

680
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
      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);
        }
700
      }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
701
      fclose(file_thread);
702
#endif  // WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
703
    }
704
#endif  // SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
705
  }
706
707

/* Print the values of the runner histogram. */
708
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
709
710
711
712
713
714
715
  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]);
716
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
717

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
718
  /* Write final output. */
719
  engine_drift_all(&e);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
720
  engine_dump_snapshot(&e);
721

722
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
723
724
  if ((res = MPI_Finalize()) != MPI_SUCCESS)
    error("call to MPI_Finalize failed with error %i.", res);
725
#endif
726

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
727
728
729
  /* Clean everything */
  engine_clean(&e);
  free(params);
730

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
731
732
  /* Say goodbye. */
  if (myrank == 0) message("done. Bye.");
733

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
734
735
  /* All is calm, all is bright. */
  return 0;
736
}