HBT.cpp

#include <cstdlib>
#include <iostream>
#include <omp.h>
#include <string>

#include "src/config_parser.h"
#include "src/datatypes.h"
#include "src/halo.h"
#include "src/mpi_wrapper.h"
#include "src/mymath.h"
#include "src/particle_exchanger.h"
#include "src/snapshot.h"
#include "src/subhalo.h"
#include "src/merger_tree.h"

#include "git_version_info.h"

int main(int argc, char **argv)
{

  MPI_Init(&argc, &argv);
  MpiWorker_t world(MPI_COMM_WORLD);

#ifdef _OPENMP
  // omp_set_nested(0);
  omp_set_max_active_levels(1); // max_active_level 0: no para; 1: single layer; >1: nest enabled
#endif

  int snapshot_start, snapshot_end;
  if (world.rank() == 0)
  {
    // Print information about the version being run.
    std::cout << "HBT-HERONS compiled using git branch: " << branch_name << " and commit: " << commit_hash;
    if (uncommitted_changes)
      std::cout << " (with uncommitted changes)";
    else
      std::cout << " (clean)";
    std::cout << std::endl;

    ParseHBTParams(argc, argv, HBTConfig, snapshot_start, snapshot_end);
    mkdir(HBTConfig.SubhaloPath.c_str(), 0755);

    std::cout << argv[0] << " run using " << world.size() << " mpi tasks";
#ifdef _OPENMP
#pragma omp parallel
#pragma omp master
    std::cout << ", each with " << omp_get_num_threads() << " threads";
    std::cout << std::endl;
#endif
    std::cout << std::endl;
    std::cout << "Configured with the following data type sizes (bytes):" << std::endl;
    std::cout << "  Real quantities    : " << sizeof(HBTReal) << std::endl;
    std::cout << "  Integer quantities : " << sizeof(HBTInt) << std::endl;
    std::cout << "  Particle velocities: " << sizeof(HBTVelType) << std::endl;
    std::cout << "  Particle masses    : " << sizeof(HBTMassType) << std::endl;
    std::cout << "  Size of Particle_t : " << sizeof(Particle_t) << std::endl;
  }
  HBTConfig.BroadCast(world, 0, snapshot_start, snapshot_end);

  SubhaloSnapshot_t subsnap;
  subsnap.Load(world, snapshot_start - 1, SubReaderDepth_t::SrcParticles);

  if (world.rank() == 0)
    std::cout << std::endl;

  /* Main loop, iterate over chosen data outputs */
  for (int isnap = snapshot_start; isnap <= snapshot_end; isnap++)
  {
    /* We start the timer for the current output, and store the time to use as a
     * reference if we want to compute fine-grained timings. */
    global_timer.Tick("start", world.Communicator);
    ReferenceTime() = global_timer.tickers[0];

    /* Load particle information */
    ParticleSnapshot_t partsnap(world, isnap);

    /* Load FOF group information. We pass partsnap because some simulation
     * formats use snapshot particles to create halo objects (e.g. GADGET4). */
    HaloSnapshot_t halosnap(world, partsnap);
    global_timer.Tick("read_halo", world.Communicator);

    /* For some input formats we load parameters directly from the snapshots,
     * so we delay writing Parameters.log until the values are known. */
    if ((isnap == snapshot_start) && (world.rank() == 0))
      HBTConfig.DumpParameters();

    halosnap.UpdateParticles(world, partsnap);
    global_timer.Tick("update_halo", world.Communicator);

    subsnap.SetSnapshotIndex(isnap);
    subsnap.UpdateParticles(world, partsnap);
    subsnap.UpdateMostBoundPosition(world, partsnap);
    global_timer.Tick("update_subhalo", world.Communicator);

    // Don't need the particle data after this point, so save memory
    partsnap.ClearParticles();

    /* Clean up the source subhaloes from duplicate particles originating from the
     * previous snapshot. We need to do it here so that any removed bound particles
     * contribute to the estimate of the subgroup CoM position and velocity (used in
     * decide centrals). We do it before assign hosts since subhaloes can change FOF
     * and ranks, making the masking difficult. */
    subsnap.CleanTracks();
    global_timer.Tick("clean_tracks", world.Communicator);

    /* We assign a FOF host to every pre-existing subhalo. All particles belonging to a
     * secondary subhalo are constrained to be within the FOF assigned to the
     * subhalo they belong to. Constraint not applied if particles are fof-less.*/
    subsnap.AssignHosts(world, halosnap, partsnap);
    global_timer.Tick("assign_hosts", world.Communicator);

    /* Store the NumTracersForDescendants most bound particles of subhaloes
     * resolved in the previous output. These will be used after unbinding to
     * determine which subhalo has accreted them. Need to do here since
     * AssignHosts will mask out some particles, and hence change the Particle
     * vector of subhaloes. Currently only correct if NumTracersForDescendants
     * <= MinNumTracerPartOfSub, because CleanTracks may have removed particles. */
    MergerTreeInfo merger_tree;
    merger_tree.StoreTracerIds(subsnap.Subhalos, HBTConfig.NumTracersForDescendants);
    global_timer.Tick("store_tracers", world.Communicator);

    /* We constrain particles to belong to FOF that hosts the subhalo they are
     * associated to. Need to do after StoringTracerIds, since this step can lead
     * to the loss of some of the most bound tracer particles */
    subsnap.ConstrainToSingleHost(halosnap);

    /* We decide which subhaloes are the central of each FOF group. Centrals are
     * assigned all the particles in the FOF that do not belong to secondary
     * subhaloes. */
    subsnap.PrepareCentrals(world, halosnap);
    global_timer.Tick("prepare_centrals", world.Communicator);

    /* Assign gas particles to the same subhalo as their nearest neighbour
       tracer type particle in the same FoF group */
    if (world.rank() == 0)
      std::cout << "Reassigning particles...\n";
    subsnap.ReassignParticles(world, halosnap);
    global_timer.Tick("reassign_particles", world.Communicator);

    /* We recursively unbind subhaloes in a depth-first approach, defined
     * by hierarchical relationships. After unbinding a given object, we check
     * wheteher any of its deeper subhaloes overlap in phase-space (if so, this
     * triggers re-unbinding). We also truncate the source of each
     * subhalo based on its number of bound particles.  */
    if (world.rank() == 0)
      std::cout << "Unbinding... ";
    subsnap.RefineParticles(world);
    global_timer.Tick("unbind", world.Communicator);

    /* Assign a unique TrackId to newly created subgroups. Update depth values,
     * hierarchical relationship, globalise FOF host values and compute other
     * subhalo properties (e.g. Vmax) */
    subsnap.UpdateTracks(world, halosnap);
    global_timer.Tick("update_tracks", world.Communicator);

    /* We locate where the tagged particles of previously bound subhaloes have
     * ended up in. */
    merger_tree.FindDescendants(subsnap.Subhalos, world);
    global_timer.Tick("merger_tree", world.Communicator);

    /* Save subhaloes */
    subsnap.Save(world);
    global_timer.Tick("write_subhalos", world.Communicator);

    /* Save timing information with root rank and reset timer for all ranks. */
    if (world.rank() == 0)
    {
      double TotalTime = global_timer.SaveSnapshotTiming(HBTConfig.SubhaloPath, subsnap.GetSnapshotIndex(), subsnap.GetSnapshotId());

      std::cout << "Snapshot " << subsnap.GetSnapshotId() << " (SnapshotIndex = " << isnap << ")" << " done. It took " << TotalTime << " seconds." << std::endl;
      std::cout << std::endl;
    }
    global_timer.Reset();
  }

  MPI_Finalize();
  return 0;
}