main.c 24.5 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
  /* How large are the parts? */
  if (myrank == 0) {
314
315
316
317
318
319
320
    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(task)      is %4zi bytes.", sizeof(struct task));
    message("sizeof(cell)      is %4zi bytes.", sizeof(struct cell));
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
321
  }
322

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
323
324
325
326
327
328
329
330
331
  /* 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");
  }
332
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
333
334
  /* Broadcast the parameter file */
  MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
335
#endif
336

337
/* Prepare the domain decomposition scheme */
338
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
339
  struct partition initial_partition;
340
  enum repartition_type reparttype;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
341
342
343
344
345
346
347
348
349
350
351
  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]);
  }
352
#endif
353

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
354
  /* Initialize unit system and constants */
355
  struct unit_system us;
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
356
357
358
359
360
361
362
363
364
365
366
367
368
369
  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
370
  if (with_hydro) hydro_props_init(&hydro_properties, params);
371

372
373
374
375
  /* 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
376
377
378
  /* Read particles and space information from (GADGET) ICs */
  char ICfileName[200] = "";
  parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
379
380
  const int replicate =
      parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
381
382
383
  if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
  fflush(stdout);

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

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

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

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

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

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
480
481
482
483
484
485
  /* 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]);
  }
486

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

  /* 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,
518
519
              engine_policies, talking, &us, &prog_const, &hydro_properties,
              &gravity_properties, &potential, &cooling_func, &sourceterms);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
520
521
522
523
524
525
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
526

527
528
/* Init the runner history. */
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
529
  for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
530
531
#endif

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

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

#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
573
574
575
  /* Split the space. */
  engine_split(&e, &initial_partition);
  engine_redistribute(&e);
576
577
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
578
579
  /* Initialise the particles */
  engine_init_particles(&e, flag_entropy_ICs);
580

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

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

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

593
594
595
596
597
/* Repartition the space amongst the nodes? */
#ifdef WITH_MPI
    if (j % 100 == 2) e.forcerepart = reparttype;
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
598
599
    /* Reset timers */
    timers_reset(timers_mask_all);
600

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
601
602
    /* Take a step. */
    engine_step(&e);
603

604
#ifdef SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
605
606
    /* Dump the task data using the given frequency. */
    if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
607
608
#ifdef WITH_MPI

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
609
610
611
612
613
614
      /* 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");
615
        fclose(file_thread);
616
      }
617
      MPI_Barrier(MPI_COMM_WORLD);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
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
654
655
656
657
658

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

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

/* Print the values of the runner histogram. */
688
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
689
690
691
692
693
694
695
  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]);
696
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
697

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
698
699
  /* Write final output. */
  engine_dump_snapshot(&e);
700

701
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
702
703
  if ((res = MPI_Finalize()) != MPI_SUCCESS)
    error("call to MPI_Finalize failed with error %i.", res);
704
#endif
705

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
706
707
708
  /* Clean everything */
  engine_clean(&e);
  free(params);
709

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
710
711
  /* Say goodbye. */
  if (myrank == 0) message("done. Bye.");
712

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
713
714
  /* All is calm, all is bright. */
  return 0;
715
}