main_fof.c 25.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
 *                    Angus Lepper (angus.lepper@ed.ac.uk)
 *               2016 John A. Regan (john.a.regan@durham.ac.uk)
 *                    Tom Theuns (tom.theuns@durham.ac.uk)
 *
 * 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.
 *
 * 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.
 *
 * 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/>.
 *
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
#include <errno.h>
#include <fenv.h>
#include <libgen.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
#include <unistd.h>

/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
#endif

/* Local headers. */
44
#include "argparse.h"
45
46
47
48
49
50
51
52
53
54
#include "swift.h"

/* Engine policy flags. */
#ifndef ENGINE_POLICY
#define ENGINE_POLICY engine_policy_none
#endif

/* Global profiler. */
struct profiler prof;

55
/*  Usage string. */
56
57
58
59
60
static const char *const fof_usage[] = {
    "fof [options] [[--] param-file]",
    "fof [options] param-file",
    "fof_mpi [options] [[--] param-file]",
    "fof_mpi [options] param-file",
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
    NULL,
};

/* Function to handle multiple -P arguments. */
struct cmdparams {
  const char *param[PARSER_MAX_NO_OF_PARAMS];
  int nparam;
};

static int handle_cmdparam(struct argparse *self,
                           const struct argparse_option *opt) {
  struct cmdparams *cmdps = (struct cmdparams *)opt->data;
  cmdps->param[cmdps->nparam] = *(char **)opt->value;
  cmdps->nparam++;
  return 1;
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
}

/**
 * @brief Main routine that loads a few particles and generates some output.
 *
 */
int main(int argc, char *argv[]) {

  struct clocks_time tic, toc;
  struct engine e;

  /* Structs used by the engine. Declare now to make sure these are always in
   * scope.  */
  struct cosmology cosmo;
  struct pm_mesh mesh;
  struct gpart *gparts = NULL;
  struct gravity_props gravity_properties;
93
  struct fof_props fof_properties;
94
95
96
97
  struct part *parts = NULL;
  struct phys_const prog_const;
  struct space s;
  struct spart *sparts = NULL;
98
  struct bpart *bparts = NULL;
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
  struct unit_system us;

  int nr_nodes = 1, myrank = 0;

#ifdef WITH_MPI
  /* Start by initializing MPI. */
  int res = 0, prov = 0;
  if ((res = MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &prov)) !=
      MPI_SUCCESS)
    error("Call to MPI_Init failed with error %i.", res);
  if (prov != MPI_THREAD_MULTIPLE)
    error(
        "MPI does not provide the level of threading"
        " required (MPI_THREAD_MULTIPLE).");
  if ((res = MPI_Comm_size(MPI_COMM_WORLD, &nr_nodes)) != MPI_SUCCESS)
    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);

118
  /* Make sure messages are stamped with the correct rank and step. */
119
  engine_rank = myrank;
120
  engine_current_step = 0;
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136

  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);
  if (myrank == 0)
    printf("[0000] [00000.0] main: MPI is up and running with %i node(s).\n\n",
           nr_nodes);
  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.");
  }
  fflush(stdout);

#endif

  /* Welcome to SWIFT, you made the right choice */
137
  if (myrank == 0) greetings(/*fof=*/1);
138
139
140
141
142

  int with_aff = 0;
  int dump_tasks = 0;
  int dump_threadpool = 0;
  int with_fp_exceptions = 0;
143
  int with_cosmology = 0;
144
145
146
  int with_stars = 0;
  int with_black_holes = 0;
  int with_hydro = 0;
147
148
  int verbose = 0;
  int nr_threads = 1;
149
150
151
  char *output_parameters_filename = NULL;
  char *cpufreqarg = NULL;
  char *param_filename = NULL;
152
  unsigned long long cpufreq = 0;
153
154
155
156
157
158
159
160
161
  struct cmdparams cmdps;
  cmdps.nparam = 0;
  cmdps.param[0] = NULL;
  char *buffer = NULL;

  /* Parse the command-line parameters. */
  struct argparse_option options[] = {
      OPT_HELP(),

162
      OPT_GROUP("  Friends-of-Friends options:\n"),
163
164
      OPT_BOOLEAN('c', "cosmology", &with_cosmology,
                  "Run with cosmological information.", NULL, 0, 0),
165
166
167
      OPT_BOOLEAN(0, "hydro", &with_hydro, "Read gas particles from the ICs.",
                  NULL, 0, 0),
      OPT_BOOLEAN(0, "stars", &with_stars, "Read stars from the ICs.", NULL, 0,
168
                  0),
169
170
      OPT_BOOLEAN(0, "black-holes", &with_black_holes,
                  "Read black holes from the ICs.", NULL, 0, 0),
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201

      OPT_GROUP("  Control options:\n"),
      OPT_BOOLEAN('a', "pin", &with_aff,
                  "Pin runners using processor affinity.", NULL, 0, 0),
      OPT_BOOLEAN('e', "fpe", &with_fp_exceptions,
                  "Enable floating-point exceptions (debugging mode).", NULL, 0,
                  0),
      OPT_STRING('f', "cpu-frequency", &cpufreqarg,
                 "Overwrite the CPU "
                 "frequency (Hz) to be used for time measurements.",
                 NULL, 0, 0),
      OPT_STRING('P', "param", &buffer,
                 "Set parameter value, overiding the value read from the "
                 "parameter file. Can be used more than once {sec:par:value}.",
                 handle_cmdparam, (intptr_t)&cmdps, 0),
      OPT_INTEGER('t', "threads", &nr_threads,
                  "The number of threads to use on each MPI rank. Defaults to "
                  "1 if not specified.",
                  NULL, 0, 0),
      OPT_INTEGER('v', "verbose", &verbose,
                  "Run in verbose mode, in MPI mode 2 outputs from all ranks.",
                  NULL, 0, 0),
      OPT_INTEGER('y', "task-dumps", &dump_tasks,
                  "Time-step frequency at which task graphs are dumped.", NULL,
                  0, 0),
      OPT_INTEGER('Y', "threadpool-dumps", &dump_threadpool,
                  "Time-step frequency at which threadpool tasks are dumped.",
                  NULL, 0, 0),
      OPT_END(),
  };
  struct argparse argparse;
202
  argparse_init(&argparse, options, fof_usage, 0);
203
204
205
206
207
208
209
210
211
212
213
214
  argparse_describe(&argparse, "\nParameters:",
                    "\nSee the file examples/parameter_example.yml for an "
                    "example of parameter file.");
  int nargs = argparse_parse(&argparse, argc, (const char **)argv);

  /* Need a parameter file. */
  if (nargs != 1) {
    if (myrank == 0) argparse_usage(&argparse);
    printf("\nError: no parameter file was supplied.\n");
    return 1;
  }
  param_filename = argv[0];
215

216
217
218
219
220
221
  /* Checks of options. */
#if !defined(HAVE_SETAFFINITY) || !defined(HAVE_LIBNUMA)
  if (with_aff) {
    printf("Error: no NUMA support for thread affinity\n");
    return 1;
  }
222
#endif
223
224
225
226
227
228

#ifndef HAVE_FE_ENABLE_EXCEPT
  if (with_fp_exceptions) {
    printf("Error: no support for floating point exceptions\n");
    return 1;
  }
229
#endif
230

231
#ifndef SWIFT_DEBUG_TASKS
232
233
234
235
236
237
238
239
  if (dump_tasks) {
    if (myrank == 0) {
      message(
          "WARNING: complete task dumps are only created when "
          "configured with --enable-task-debugging.");
      message("         Basic task statistics will be output.");
    }
  }
240
#endif
241

242
#ifndef SWIFT_DEBUG_THREADPOOL
243
244
245
246
247
248
  if (dump_threadpool) {
    printf(
        "Error: threadpool dumping is only possible if SWIFT was "
        "configured with the --enable-threadpool-debugging option.\n");
    return 1;
  }
249
#endif
250
251
252
253
254
255
256

  /* The CPU frequency is a long long, so we need to parse that ourselves. */
  if (cpufreqarg != NULL) {
    if (sscanf(cpufreqarg, "%llu", &cpufreq) != 1) {
      if (myrank == 0)
        printf("Error parsing CPU frequency (%s).\n", cpufreqarg);
      return 1;
257
    }
258
  }
259
260

  /* Write output parameter file */
261
  if (myrank == 0 && output_parameters_filename != NULL) {
262
    io_write_output_field_parameter(output_parameters_filename, with_cosmology);
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
    printf("End of run.\n");
    return 0;
  }

/* Let's pin the main thread, now we know if affinity will be used. */
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
  if (with_aff &&
      ((ENGINE_POLICY)&engine_policy_setaffinity) == engine_policy_setaffinity)
    engine_pin();
#endif

  /* Genesis 1.1: And then, there was time ! */
  clocks_set_cpufreq(cpufreq);

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

  /* Report CPU frequency.*/
  cpufreq = clocks_get_cpufreq();
  if (myrank == 0) {
    message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
  }

/* Report host name(s). */
#ifdef WITH_MPI
  if (talking) {
    message("Rank %d running on: %s", myrank, hostname());
  }
#else
  message("Running on: %s", hostname());
#endif

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

301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
  /* Do we choke on FP-exceptions ? */
  if (with_fp_exceptions) {
#ifdef HAVE_FE_ENABLE_EXCEPT
    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
    if (myrank == 0)
      message("WARNING: Floating point exceptions will be reported.");
  }

/* Do we have slow barriers? */
#ifndef HAVE_PTHREAD_BARRIERS
  if (myrank == 0)
    message("WARNING: Non-optimal thread barriers are being used.");
#endif

  /* How large are the parts? */
  if (myrank == 0) {
    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));
  }

  /* Read the parameter file */
  struct swift_params *params =
      (struct swift_params *)malloc(sizeof(struct swift_params));
  if (params == NULL) error("Error allocating memory for the parameter file.");
  if (myrank == 0) {
333
334
    message("Reading runtime parameters from file '%s'", param_filename);
    parser_read_file(param_filename, params);
335
336

    /* Handle any command-line overrides. */
337
    if (cmdps.nparam > 0) {
338
339
340
      message(
          "Overwriting values read from the YAML file with command-line "
          "values.");
341
342
      for (int k = 0; k < cmdps.nparam; k++)
        parser_set_param(params, cmdps.param[k]);
343
344
345
346
347
348
    }
  }
#ifdef WITH_MPI
  /* Broadcast the parameter file */
  MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
#endif
349

350
351
352
353
354
355
356
  /* "Read" the Select Output file - this should actually do nothing, but
   * we need to mock the struct up for passing to `engine_init` */

  struct output_options *output_options =
      (struct output_options *)malloc(sizeof(struct output_options));
  output_options_init(params, myrank, output_options);

357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
  /* Check that we can write the snapshots by testing if the output
   * directory exists and is searchable and writable. */
  char basename[PARSER_MAX_LINE_SIZE];
  parser_get_param_string(params, "Snapshots:basename", basename);
  const char *dirp = dirname(basename);
  if (access(dirp, W_OK | X_OK) != 0) {
    error("Cannot write snapshots in directory %s (%s)", dirp, strerror(errno));
  }

  /* Prepare the domain decomposition scheme */
  struct repartition reparttype;
#ifdef WITH_MPI
  struct partition initial_partition;
  partition_init(&initial_partition, &reparttype, params, nr_nodes);

  /* Let's report what we did */
  if (myrank == 0) {
#if defined(HAVE_PARMETIS)
    if (reparttype.usemetis)
      message("Using METIS serial partitioning:");
    else
      message("Using ParMETIS partitioning:");
379
#elif defined(HAVE_METIS)
380
    message("Using METIS serial partitioning:");
381
382
#else
    message("Non-METIS partitioning:");
383
384
385
386
387
388
389
390
391
392
#endif
    message("  initial partitioning: %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("  repartitioning: %s", repartition_name[reparttype.type]);
  }
#endif

393
394
395
396
  /* Prepare and verify the selection of outputs */
  io_prepare_output_fields(output_options, with_cosmology, /*with_fof=*/1,
                           /*with_structure_finding=*/0);

397
398
399
  /* Initialize unit system and constants */
  units_init_from_params(&us, params, "InternalUnitSystem");
  phys_const_init(&us, params, &prog_const);
400
  if (myrank == 0) {
401
402
403
404
405
406
    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);
407
  }
408

409
410
411
412
  /* Read particles and space information from ICs */
  char ICfileName[200] = "";
  parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
  const int periodic =
413
      parser_get_param_int(params, "InitialConditions:periodic");
414
  const int replicate =
415
      parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
416
  const int clean_smoothing_length_values = parser_get_opt_param_int(
417
      params, "InitialConditions:cleanup_smoothing_lengths", 0);
418
  const int cleanup_h = parser_get_opt_param_int(
419
      params, "InitialConditions:cleanup_h_factors", 0);
420
  const int cleanup_sqrt_a = parser_get_opt_param_int(
421
      params, "InitialConditions:cleanup_velocity_factors", 0);
422

423
  /* Initialise the cosmology */
424
425
426
427
428
429
430
431
432
433
434
  if (with_cosmology)
    cosmology_init(params, &us, &prog_const, &cosmo);
  else
    cosmology_init_no_cosmo(&cosmo);
  if (myrank == 0 && with_cosmology) cosmology_print(&cosmo);

  /* Initialise the equation of state */
  if (with_hydro)
    eos_init(&eos, &prog_const, &us, params);
  else
    bzero(&eos, sizeof(struct eos_parameters));
435

436
437
  /* Initialise the FOF properties */
  bzero(&fof_properties, sizeof(struct fof_props));
438
  fof_init(&fof_properties, params, &prog_const, &us, /*stand-alone=*/1);
439

440
441
442
443
444
445
446
  /* Be verbose about what happens next */
  if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
  if (myrank == 0 && cleanup_h)
    message("Cleaning up h-factors (h=%f)", cosmo.h);
  if (myrank == 0 && cleanup_sqrt_a)
    message("Cleaning up a-factors from velocity (a=%f)", cosmo.a);
  fflush(stdout);
447

448
  /* Get ready to read particles of all kinds */
449
  int flag_entropy_ICs = 0;
450
  size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
451
  double dim[3] = {0., 0., 0.};
452

453
  if (myrank == 0) clocks_gettime(&tic);
454
455
456
#if defined(HAVE_HDF5)
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
457
458
459
460
461
462
463
  read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
                   &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
                   &flag_entropy_ICs, with_hydro,
                   /*with_grav=*/1, with_stars, with_black_holes,
                   with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
                   myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
                   /*dry_run=*/0);
464
#else
465
466
467
468
469
470
  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas,
                 &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
                 &flag_entropy_ICs, with_hydro,
                 /*with_grav=*/1, with_stars, with_black_holes, with_cosmology,
                 cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes,
                 MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, /*dry_run=*/0);
471
472
#endif
#else
Josh Borrow's avatar
Josh Borrow committed
473
474
475
  read_ic_single(
      ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart,
      &Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
476
477
      /*with_grav=*/1, with_stars, with_black_holes, with_cosmology, cleanup_h,
      cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0);
478
479
#endif
#endif
480
481
  if (myrank == 0) {
    clocks_gettime(&toc);
482
483
    message("Reading initial conditions took %.3f %s.", clocks_diff(&tic, &toc),
            clocks_getunit());
484
485
    fflush(stdout);
  }
486

487
#ifdef SWIFT_DEBUG_CHECKS
488
  /* Check that the other links are correctly set */
489
490
  part_verify_links(parts, gparts, sparts, bparts, Ngas, Ngpart, Nspart, Nbpart,
                    1);
491
#endif
492

493
  /* Get the total number of particles across all nodes. */
494
  long long N_total[swift_type_count + 1] = {0};
495
  long long Nbaryons = Ngas + Nspart + Nbpart;
496
#if defined(WITH_MPI)
497
498
499
500
501
502
503
504
505
  long long N_long[swift_type_count + 1] = {0};
  N_long[swift_type_gas] = Ngas;
  N_long[swift_type_dark_matter] = Ngpart - Ngpart_background - Nbaryons;
  N_long[swift_type_dark_matter_background] = Ngpart_background;
  N_long[swift_type_stars] = Nspart;
  N_long[swift_type_black_hole] = Nbpart;
  N_long[swift_type_count] = Ngpart;
  MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
                MPI_SUM, MPI_COMM_WORLD);
506
#else
507
508
509
510
511
512
  N_total[swift_type_gas] = Ngas;
  N_total[swift_type_dark_matter] = Ngpart - Ngpart_background - Nbaryons;
  N_total[swift_type_dark_matter_background] = Ngpart_background;
  N_total[swift_type_stars] = Nspart;
  N_total[swift_type_black_hole] = Nbpart;
  N_total[swift_type_count] = Ngpart;
513
#endif
514

515
516
  if (myrank == 0)
    message(
517
        "Read %lld gas particles, %lld stars particles, %lld black hole "
518
519
520
521
522
523
        "particles, %lld DM particles and %lld DM background particles from "
        "the ICs.",
        N_total[swift_type_gas], N_total[swift_type_stars],
        N_total[swift_type_black_hole], N_total[swift_type_dark_matter],
        N_total[swift_type_dark_matter_background]);

524
525
526
527
528
  const int with_DM_particles = N_total[swift_type_dark_matter] > 0;
  const int with_baryon_particles =
      (N_total[swift_type_gas] + N_total[swift_type_stars] +
       N_total[swift_type_black_hole]) > 0;

529
530
531
  /* Do we have background DM particles? */
  const int with_DM_background_particles =
      N_total[swift_type_dark_matter_background] > 0;
532

533
534
  /* Initialize the space with these data. */
  if (myrank == 0) clocks_gettime(&tic);
535
536
  space_init(&s, params, &cosmo, dim, /*hydro_props=*/NULL, parts, gparts,
             sparts, bparts, Ngas, Ngpart, Nspart, Nbpart, periodic, replicate,
537
             /*generate_gas_in_ics=*/0, /*hydro=*/N_total[0] > 0, /*gravity=*/1,
538
             /*with_star_formation=*/0, with_DM_background_particles, talking,
539
             /*dry_run=*/0, nr_nodes);
540

541
542
543
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
544
            clocks_getunit());
545
546
    fflush(stdout);
  }
547

548
549
550
  /* Initialise the gravity scheme */
  bzero(&gravity_properties, sizeof(struct gravity_props));
  gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
551
                     with_cosmology, /*with_external_gravity=*/0,
552
                     with_baryon_particles, with_DM_particles,
553
                     with_DM_background_particles, periodic, s.dim);
554

555
556
  /* Initialise the long-range gravity mesh */
  if (periodic) {
557
#ifdef HAVE_FFTW
558
    pm_mesh_init(&mesh, &gravity_properties, s.dim, nr_threads);
559
#else
560
    /* Need the FFTW library if periodic and self gravity. */
561
    error("No FFTW library found. Cannot compute periodic long-range forces.");
562
#endif
563
564
565
  } else {
    pm_mesh_init_no_mesh(&mesh, s.dim);
  }
566

567
568
  /* Also update the total counts (in case of changes due to replication) */
  Nbaryons = s.nr_parts + s.nr_sparts + s.nr_bparts;
569
#if defined(WITH_MPI)
570
571
572
573
574
575
576
  N_long[swift_type_gas] = s.nr_parts;
  N_long[swift_type_dark_matter] = s.nr_gparts - Ngpart_background - Nbaryons;
  N_long[swift_type_count] = s.nr_gparts;
  N_long[swift_type_stars] = s.nr_sparts;
  N_long[swift_type_black_hole] = s.nr_bparts;
  MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
                MPI_SUM, MPI_COMM_WORLD);
577
#else
578
579
580
581
582
  N_total[swift_type_gas] = s.nr_parts;
  N_total[swift_type_dark_matter] = s.nr_gparts - Ngpart_background - Nbaryons;
  N_total[swift_type_count] = s.nr_gparts;
  N_total[swift_type_stars] = s.nr_sparts;
  N_total[swift_type_black_hole] = s.nr_bparts;
583
#endif
584

585
586
587
  /* 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],
588
            s.dim[2]);
589
590
    message("space %s periodic.", s.periodic ? "is" : "isn't");
    message("highest-level cell dimensions are [ %i %i %i ].", s.cdim[0],
591
            s.cdim[1], s.cdim[2]);
592
593
594
595
596
    message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
    message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
    message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
    message("maximum depth is %d.", s.maxdepth);
    fflush(stdout);
597
  }
598

599
600
601
602
  /* Construct the engine policy */
  int engine_policies = ENGINE_POLICY | engine_policy_steal;
  engine_policies |= engine_policy_self_gravity;
  engine_policies |= engine_policy_fof;
603
  if (with_cosmology) engine_policies |= engine_policy_cosmology;
604

605
606
  /* Initialize the engine with the space and policies. */
  if (myrank == 0) clocks_gettime(&tic);
607
608
609
610
611
612
613
614
615
616
617
  engine_init(&e, &s, params, output_options, N_total[swift_type_gas],
              N_total[swift_type_count], N_total[swift_type_stars],
              N_total[swift_type_black_hole],
              N_total[swift_type_dark_matter_background], engine_policies,
              talking, &reparttype, &us, &prog_const, &cosmo,
              /*hydro_properties=*/NULL, /*entropy_floor=*/NULL,
              &gravity_properties,
              /*stars_properties=*/NULL, /*black_holes_properties=*/NULL,
              /*feedback_properties=*/NULL, &mesh, /*potential=*/NULL,
              /*cooling_func=*/NULL, /*starform=*/NULL, /*chemistry=*/NULL,
              &fof_properties, /*los_properties=*/NULL);
618
619
620
  engine_config(/*restart=*/0, /*fof=*/1, &e, params, nr_nodes, myrank,
                nr_threads, with_aff, talking, NULL);

621
622
623
  if (myrank == 0) {
    clocks_gettime(&toc);
    message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
624
            clocks_getunit());
625
626
    fflush(stdout);
  }
627

628
629
  /* Get some info to the user. */
  if (myrank == 0) {
630
631
    const long long N_DM = N_total[swift_type_dark_matter] +
                           N_total[swift_type_dark_matter_background];
632
    message(
633
        "Running FOF on %lld gas particles, %lld stars particles %lld black "
634
        "hole particles and %lld DM particles (%lld gravity particles)",
635
636
        N_total[swift_type_gas], N_total[swift_type_stars],
        N_total[swift_type_black_hole], N_DM, N_total[swift_type_count]);
637
    message(
638
639
640
641
        "from t=%.3e until t=%.3e with %d ranks, %d threads / rank and %d "
        "task queues / rank (dt_min=%.3e, dt_max=%.3e)...",
        e.time_begin, e.time_end, nr_nodes, e.nr_threads, e.sched.nr_queues,
        e.dt_min, e.dt_max);
642
643
    fflush(stdout);
  }
644

645
#ifdef WITH_MPI
646
  /* Split the space. */
647
648
  engine_split(&e, &initial_partition);
  engine_redistribute(&e);
649
#endif
650

James Willis's avatar
James Willis committed
651
#ifdef SWIFT_DEBUG_TASKS
652
  e.tic_step = getticks();
James Willis's avatar
James Willis committed
653
#endif
654

655
656
  /* Initialise the tree and communication tasks */
  engine_rebuild(&e, /*repartitionned=*/0, clean_smoothing_length_values);
657

James Willis's avatar
James Willis committed
658
#ifdef SWIFT_DEBUG_TASKS
659
  e.toc_step = getticks();
James Willis's avatar
James Willis committed
660
#endif
661

662
663
  /* Perform the FOF search */
  engine_fof(&e, /*dump_results=*/1, /*seed_black_holes=*/0);
James Willis's avatar
James Willis committed
664

665
666
  /* Write output. */
  engine_dump_snapshot(&e);
667

668
#ifdef WITH_MPI
669
  MPI_Barrier(MPI_COMM_WORLD);
670
#endif
671

672
673
  /* Dump the task data using the given frequency. */
  if (dump_tasks) {
James Willis's avatar
James Willis committed
674
#ifdef SWIFT_DEBUG_TASKS
675
676
    task_dump_all(&e, 0);
#endif
677

678
679
680
    /* Generate the task statistics. */
    char dumpfile[40];
    snprintf(dumpfile, 40, "thread_stats-step%d.dat", 0);
681
682
    task_dump_stats(dumpfile, &e, /* dump_tasks_threshold = */ 0.f,
                    /* header = */ 0, /* allranks = */ 1);
James Willis's avatar
James Willis committed
683
  }
684

685
686
687
#ifdef SWIFT_DEBUG_THREADPOOL
  /* Dump the task data using the given frequency. */
  if (dump_threadpool) {
James Willis's avatar
James Willis committed
688
    char threadpool_dumpfile[52];
689
#ifdef WITH_MPI
690
691
    snprintf(threadpool_dumpfile, 52, "threadpool_info-rank%d-step%d.dat",
             engine_rank, 0);
692
#else
James Willis's avatar
James Willis committed
693
    snprintf(threadpool_dumpfile, 52, "threadpool_info-step%d.dat", 0);
694
#endif  // WITH_MPI
James Willis's avatar
James Willis committed
695
    threadpool_dump_log(&e.threadpool, threadpool_dumpfile, 1);
696
697
698
699
700
  } else {
    threadpool_reset_log(&e.threadpool);
  }
#endif  // SWIFT_DEBUG_THREADPOOL

701
  /* used parameters */
702
  parser_write_params_to_file(params, "fof_used_parameters.yml", /*used=*/1);
703
  /* unused parameters */
704
  parser_write_params_to_file(params, "fof_unused_parameters.yml", /*used=*/0);
705

706
707
708
709
710
711
712
713
714
715
716
717
718
  /* Dump memory use report */
#ifdef SWIFT_MEMUSE_REPORTS
  {
    char dumpfile[40];
#ifdef WITH_MPI
    snprintf(dumpfile, 40, "memuse_report-rank%d-fof.dat", engine_rank);
#else
    snprintf(dumpfile, 40, "memuse_report-fof.dat");
#endif  // WITH_MPI
    memuse_log_dump(dumpfile);
  }
#endif

719
720
721
722
723
#ifdef WITH_MPI
  if ((res = MPI_Finalize()) != MPI_SUCCESS)
    error("call to MPI_Finalize failed with error %i.", res);
#endif

724
  /* Clean everything */
725
726
  cosmology_clean(&cosmo);
  pm_mesh_clean(&mesh);
Matthieu Schaller's avatar
Matthieu Schaller committed
727
  engine_clean(&e, /*fof=*/1, /*restart=*/0);
728
  free(params);
729
  free(output_options);
730
731
732
733
734
735
736

  /* Say goodbye. */
  if (myrank == 0) message("done. Bye.");

  /* All is calm, all is bright. */
  return 0;
}