main.c 25.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
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", "-T", "", "Print timers every time-step.");
94
  printf("  %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity.");
Matthieu Schaller's avatar
Matthieu Schaller committed
95
  printf("  %2s %8s %s\n", "", "", "1: MPI-rank 0 writes ");
96
  printf("  %2s %8s %s\n", "", "", "2: All MPI-ranks write");
97
  printf("  %2s %8s %s\n", "-y", "{int}",
98
99
         "Time-step frequency at which task graphs are dumped.");
  printf("  %2s %8s %s\n", "-h", "", "Print this help message and exit.");
100
101
102
  printf(
      "\nSee the file parameter_example.yml for an example of "
      "parameter file.\n");
103
104
}

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

111
  struct clocks_time tic, toc;
112

113
  int nr_nodes = 1, myrank = 0;
114

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

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

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

143
#endif
144

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

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

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

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

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

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

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

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

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

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

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
329
330
  /* How large are the parts? */
  if (myrank == 0) {
331
332
333
334
335
336
337
338
    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
339
  }
340

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

355
  /* Prepare the domain decomposition scheme */
356
  struct repartition reparttype;
357
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
358
359
360
361
362
363
364
365
366
367
  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]);
368
    message("Using %s repartitioning", repartition_name[reparttype.type]);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
369
  }
370
#endif
371

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

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

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

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

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

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

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

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
498
499
500
501
502
503
  /* 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]);
  }
504

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

  /* Initialize the engine with the space and policies. */
  if (myrank == 0) clocks_gettime(&tic);
  struct engine e;
535
  engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, N_total[0],
Peter W. Draper's avatar
Peter W. Draper committed
536
537
538
              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
539
540
541
542
543
544
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
    fflush(stdout);
  }
545

546
547
/* Init the runner history. */
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
548
  for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
549
550
#endif

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

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

#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
592
593
594
  /* Split the space. */
  engine_split(&e, &initial_partition);
  engine_redistribute(&e);
595
596
#endif

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
597
598
  /* Initialise the particles */
  engine_init_particles(&e, flag_entropy_ICs);
599

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
603
  /* Legend */
604
  if (myrank == 0) {
605
606
607
    printf("# %6s %14s %14s %10s %10s %10s %16s [%s]\n", "Step", "Time",
           "Time-step", "Updates", "g-Updates", "s-Updates", "Wall-clock time",
           clocks_getunit());
608

609
610
611
612
613
614
615
    if (with_verbose_timers) {
      printf("timers: ");
      for (int k = 0; k < timer_count; k++) printf("%s\t", timers_names[k]);
      printf("\n");
    }
  }

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

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
619
620
    /* Reset timers */
    timers_reset(timers_mask_all);
621

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
622
623
    /* Take a step. */
    engine_step(&e);
624

625
626
627
628
629
630
631
632
633
    /* Print the timers. */
    if (with_verbose_timers) {
      printf("timers: ");
      for (int k = 0; k < timer_count; k++)
        printf("%.3f\t", clocks_from_ticks(timers[k]));
      printf("\n");
      timers_reset(0xFFFFFFFFllu);
    }

634
#ifdef SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
635
636
    /* Dump the task data using the given frequency. */
    if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
637
638
#ifdef WITH_MPI

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
639
640
641
642
643
644
      /* 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");
645
        fclose(file_thread);
646
      }
647
      MPI_Barrier(MPI_COMM_WORLD);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
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
679
680
681
682
683
684
685
686
687
688

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

690
#else
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
      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);
        }
710
      }
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
711
      fclose(file_thread);
712
#endif  // WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
713
    }
714
#endif  // SWIFT_DEBUG_TASKS
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
715
  }
716
717

/* Print the values of the runner histogram. */
718
#ifdef HIST
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
719
720
721
722
723
724
725
  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]);
726
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
727

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
728
  /* Write final output. */
729
  engine_drift_all(&e);
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
730
  engine_dump_snapshot(&e);
731

732
#ifdef WITH_MPI
Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
733
734
  if ((res = MPI_Finalize()) != MPI_SUCCESS)
    error("call to MPI_Finalize failed with error %i.", res);
735
#endif
736

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
737
738
739
  /* Clean everything */
  engine_clean(&e);
  free(params);
740

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
741
742
  /* Say goodbye. */
  if (myrank == 0) message("done. Bye.");
743

Pedro Gonnet's avatar
oops.    
Pedro Gonnet committed
744
745
  /* All is calm, all is bright. */
  return 0;
746
}