init.c

Go to the documentation of this file.

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_sf_gamma.h>
 
#include "init.h"
#include "utils.h"
 
#include "cooling.h"
#include "forcetree.h"
#include "density.h"
 
#include "timefac.h"
#include "petaio.h"
#include "domain.h"
#include "walltime.h"
#include "slotsmanager.h"
#include "hydra.h"
#include "sfr_eff.h"
#include "exchange.h"
#include "fof.h"
#include "timestep.h"
#include "timebinmgr.h"
#include "cosmology.h"
#include "gravity.h"
#include "physconst.h"
 
static struct init_params
{
    /* Gas temperature in the ICs*/
    double InitGasTemp;
    double PartAllocFactor;
} InitParams;
 
/*Set the global parameters*/
void
set_init_params(ParameterSet * ps)
{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        InitParams.InitGasTemp = param_get_double(ps, "InitGasTemp");
        InitParams.PartAllocFactor = param_get_double(ps, "PartAllocFactor");
    }
    MPI_Bcast(&InitParams, sizeof(InitParams), MPI_BYTE, 0, MPI_COMM_WORLD);
}
 
void init_timeline(int RestartSnapNum, double TimeMax, const struct header_data * header, const int SnapshotWithFOF)
{
    /*Add TimeInit and TimeMax to the output list*/
    if (RestartSnapNum < 0) {
        /* allow a first snapshot at IC time; */
        setup_sync_points(header->TimeIC, TimeMax, 0.0, SnapshotWithFOF);
    } else {
        /* skip dumping the exactly same snapshot */
        setup_sync_points(header->TimeIC, TimeMax, header->TimeSnapshot, SnapshotWithFOF);
        /* If TimeInit is not in a sensible place on the integer timeline
         * (can happen if the outputs changed since it was written)
         * start the integer timeline anew from TimeInit */
        inttime_t ti_init = ti_from_loga(log(header->TimeSnapshot)) % TIMEBASE;
        if(round_down_power_of_two(ti_init) != ti_init) {
            message(0,"Resetting integer timeline (as %x != %x) to current snapshot\n",ti_init, round_down_power_of_two(ti_init));
            setup_sync_points(header->TimeSnapshot, TimeMax, header->TimeSnapshot, SnapshotWithFOF);
        }
    }
}
 
 
static void get_mean_separation(double * MeanSeparation, const double BoxSize, const int64_t * NTotalInit);
static void check_omega(struct part_manager_type * PartManager, Cosmology * CP, int generations, double G, double * MassTable);
static void check_positions(struct part_manager_type * PartManager);
static void check_smoothing_length(struct part_manager_type * PartManager, double * MeanSpacing);
static void init_alloc_particle_slot_memory(struct part_manager_type * PartManager, struct slots_manager_type * SlotsManager, const double PartAllocFactor, struct header_data * header, MPI_Comm Comm);
 
inttime_t init(int RestartSnapNum, const char * OutputDir, struct header_data * header, Cosmology * CP)
{
    int i;
 
    init_alloc_particle_slot_memory(PartManager, SlotsManager, InitParams.PartAllocFactor, header, MPI_COMM_WORLD);
 
    /*Read the snapshot*/
    petaio_read_snapshot(RestartSnapNum, OutputDir, CP, header, PartManager, SlotsManager, MPI_COMM_WORLD);
 
    domain_test_id_uniqueness(PartManager);
 
    check_omega(PartManager, CP, get_generations(), CP->GravInternal, header->MassTable);
 
    check_positions(PartManager);
 
    double MeanSeparation[6] = {0};
 
    get_mean_separation(MeanSeparation, PartManager->BoxSize, header->NTotalInit);
 
    if(RestartSnapNum == -1)
        check_smoothing_length(PartManager, MeanSeparation);
 
    /* As the above will mostly take place
     * on Task 0, there will be a lot of imbalance*/
    MPIU_Barrier(MPI_COMM_WORLD);
 
    gravshort_set_softenings(MeanSeparation[1]);
    fof_init(MeanSeparation[1]);
 
    inttime_t Ti_Current = init_timebins(header->TimeSnapshot);
 
    #pragma omp parallel for
    for(i = 0; i < PartManager->NumPart; i++)   /* initialize sph_properties */
    {
        int j;
        P[i].Ti_drift = Ti_Current;
 
        if(RestartSnapNum == -1 && P[i].Type == 5 )
        {
            /* Note: Gadget-3 sets this to the seed black hole mass.*/
            BHP(i).Mass = P[i].Mass;
        }
 
        if(P[i].Type == 4 )
        {
            /* Touch up zero star smoothing lengths, not saved in the snapshots.*/
            if(P[i].Hsml == 0)
                P[i].Hsml = 0.1 * MeanSeparation[0];
        }
 
        if(P[i].Type == 5)
        {
            for(j = 0; j < 3; j++) {
                BHP(i).DFAccel[j] = 0;
                BHP(i).DragAccel[j] = 0;
            }
        }
 
        P[i].Key = PEANO(P[i].Pos, PartManager->BoxSize);
 
        if(P[i].Type != 0) continue;
 
        for(j = 0; j < 3; j++)
        {
            SPHP(i).HydroAccel[j] = 0;
        }
 
        if(!isfinite(SPHP(i).DelayTime ))
            endrun(6, "Bad DelayTime %g for part %d id %ld\n", SPHP(i).DelayTime, i, P[i].ID);
        SPHP(i).DtEntropy = 0;
 
        if(RestartSnapNum == -1)
        {
            SPHP(i).Density = -1;
            SPHP(i).EgyWtDensity = -1;
            SPHP(i).DhsmlEgyDensityFactor = -1;
            SPHP(i).Entropy = -1;
            SPHP(i).Ne = 1.0;
            SPHP(i).DivVel = 0;
            SPHP(i).CurlVel = 0;
            SPHP(i).DelayTime = 0;
            SPHP(i).Metallicity = 0;
            memset(SPHP(i).Metals, 0, NMETALS*sizeof(float));
            /* Initialise to primordial abundances for H and He*/
            SPHP(i).Metals[0] = HYDROGEN_MASSFRAC;
            SPHP(i).Metals[1] = 1- HYDROGEN_MASSFRAC;
            SPHP(i).Sfr = 0;
            SPHP(i).MaxSignalVel = 0;
        }
    }
    walltime_measure("/Init");
    return Ti_Current;
}
 
 
void check_omega(struct part_manager_type * PartManager, Cosmology * CP, int generations, double G, double * MassTable)
{
    double mass = 0, masstot, omega;
    int i, badmass = 0;
    int64_t totbad;
 
    #pragma omp parallel for reduction(+: mass) reduction(+: badmass)
    for(i = 0; i < PartManager->NumPart; i++) {
        /* In case zeros have been written to the saved mass array,
         * recover the true masses*/
        if(P[i].Mass == 0) {
            P[i].Mass = MassTable[P[i].Type] * ( 1. - (double)P[i].Generation/generations);
            badmass++;
        }
        mass += P[i].Mass;
    }
 
    MPI_Allreduce(&mass, &masstot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
    sumup_large_ints(1, &badmass, &totbad);
    if(totbad)
        message(0, "Warning: recovering from %ld Mass entries corrupted on disc\n",totbad);
 
    omega =
        masstot / (PartManager->BoxSize * PartManager->BoxSize * PartManager->BoxSize) / (3 * CP->Hubble * CP->Hubble / (8 * M_PI * G));
 
    /*Add the density for analytically follows massive neutrinos*/
    if(CP->MassiveNuLinRespOn)
        omega += get_omega_nu_nopart(&CP->ONu, 1);
    if(fabs(omega - CP->Omega0) > 1.0e-3)
    {
        endrun(0, "The mass content accounts only for Omega=%g,\nbut you specified Omega=%g in the parameterfile.\n",
                omega, CP->Omega0);
    }
}
 
/* Allocate the memory for particles and slots. First the total amount of particles are counted, then allocations are made*/
static void
init_alloc_particle_slot_memory(struct part_manager_type * PartManager, struct slots_manager_type * SlotsManager, const double PartAllocFactor, struct header_data * header, MPI_Comm Comm)
{
    int NTask, ThisTask;
    MPI_Comm_size(Comm, &NTask);
    MPI_Comm_rank(Comm, &ThisTask);
 
    int ptype;
 
    int64_t TotNumPart = 0;
    int64_t TotNumPartInit= 0;
    for(ptype = 0; ptype < 6; ptype ++) {
        TotNumPart += header->NTotal[ptype];
        TotNumPartInit += header->NTotalInit[ptype];
    }
 
    message(0, "Total number of particles: %018ld\n", TotNumPart);
 
    const char * PARTICLE_TYPE_NAMES [] = {"Gas", "DarkMatter", "Neutrino", "Unknown", "Star", "BlackHole"};
 
    for(ptype = 0; ptype < 6; ptype ++) {
        double MeanSeparation = header->BoxSize / pow(header->NTotalInit[ptype], 1.0 / 3);
        message(0, "% 11s: Total: %018ld Init: %018ld Mean-Sep %g \n",
                PARTICLE_TYPE_NAMES[ptype], header->NTotal[ptype], header->NTotalInit[ptype], MeanSeparation);
    }
 
    /* sets the maximum number of particles that may reside on a processor */
    int MaxPart = (int) (PartAllocFactor * TotNumPartInit / NTask);
 
    /*Allocate the particle memory*/
    particle_alloc_memory(PartManager, header->BoxSize, MaxPart);
 
    for(ptype = 0; ptype < 6; ptype ++) {
        int64_t start = ThisTask * header->NTotal[ptype] / NTask;
        int64_t end = (ThisTask + 1) * header->NTotal[ptype] / NTask;
        header->NLocal[ptype] = end - start;
        PartManager->NumPart += header->NLocal[ptype];
    }
 
    /* Allocate enough memory for stars and black holes.
     * This will be dynamically increased as needed.*/
 
    if(PartManager->NumPart >= PartManager->MaxPart) {
        endrun(1, "Overwhelmed by part: %ld > %ld\n", PartManager->NumPart, PartManager->MaxPart);
    }
 
    /* Now allocate memory for the secondary particle data arrays.
     * This may be dynamically resized later!*/
 
    /*Ensure all processors have initially the same number of particle slots*/
    int64_t newSlots[6] = {0};
 
    /* Can't use MPI_IN_PLACE, which is broken for arrays and MPI_MAX at least on intel mpi 19.0.5*/
    MPI_Allreduce(header->NLocal, newSlots, 6, MPI_INT64, MPI_MAX, Comm);
 
    for(ptype = 0; ptype < 6; ptype ++) {
            newSlots[ptype] *= PartAllocFactor;
    }
    /* Boost initial amount of stars allocated, as it is often uneven.
     * The total number of stars is usually small so this doesn't
     * waste that much memory*/
    newSlots[4] *= 2;
 
    slots_reserve(0, newSlots, SlotsManager);
 
    /* so we can set up the memory topology of secondary slots */
    slots_setup_topology(PartManager, header->NLocal, SlotsManager);
}
 
void check_positions(struct part_manager_type * PartManager)
{
    int i;
    int numzero = 0;
    int lastzero = -1;
 
    #pragma omp parallel for reduction(+: numzero) reduction(max:lastzero)
    for(i=0; i< PartManager->NumPart; i++){
        int j;
        double * Pos = PartManager->Base[i].Pos;
        for(j=0; j<3; j++) {
            if(Pos[j] < 0 || Pos[j] > PartManager->BoxSize || !isfinite(Pos[j]))
                endrun(0,"Particle %d is outside the box (L=%g) at (%g %g %g)\n",i, PartManager->BoxSize, Pos[0], Pos[1], Pos[2]);
        }
        if((Pos[0] < 1e-35) && (Pos[1] < 1e-35) && (Pos[2] < 1e-35)) {
            numzero++;
            lastzero = i;
        }
    }
    if(numzero > 1)
        endrun(5, "Particle positions contain %d zeros at particle %d. Pos %g %g %g. Likely write corruption!\n",
                numzero, lastzero, PartManager->Base[lastzero].Pos[0], PartManager->Base[lastzero].Pos[1], PartManager->Base[lastzero].Pos[2]);
}
 
void check_smoothing_length(struct part_manager_type * PartManager, double * MeanSpacing)
{
    int i;
    int numprob = 0;
    int lastprob = -1;
    #pragma omp parallel for reduction(+: numprob) reduction(max:lastprob)
    for(i=0; i< PartManager->NumPart; i++){
        if(P[i].Type != 5 && P[i].Type != 0)
            continue;
        if(P[i].Hsml > PartManager->BoxSize || P[i].Hsml <= 0) {
            P[i].Hsml = MeanSpacing[P[i].Type];
            numprob++;
            lastprob = i;
        }
    }
    if(numprob > 0)
        message(5, "Bad smoothing lengths %d last bad %d hsml %g id %ld\n", numprob, lastprob, P[lastprob].Hsml, P[lastprob].ID);
}
 
/* When we restart, validate the SPH properties of the particles.
 * This also allows us to increase MinEgySpec on a restart if we choose.*/
void check_density_entropy(Cosmology * CP, const double MinEgySpec, const double atime)
{
    const double a3 = pow(atime, 3);
    int i;
    int bad = 0;
    double meanbar = CP->OmegaBaryon * 3 * HUBBLE * CP->HubbleParam * HUBBLE * CP->HubbleParam/ (8 * M_PI * GRAVITY);
    #pragma omp parallel for reduction(+: bad)
    for(i = 0; i < SlotsManager->info[0].size; i++) {
        /* This allows us to continue gracefully if
            * there was some kind of bug in the run that output the snapshot.
            * density() below will fix this up.*/
        if(SphP[i].Density <= 0 || !isfinite(SphP[i].Density)) {
            SphP[i].Density = meanbar;
            bad++;
        }
        if(SphP[i].EgyWtDensity <= 0 || !isfinite(SphP[i].EgyWtDensity)) {
            SphP[i].EgyWtDensity = SphP[i].Density;
        }
        double minent = GAMMA_MINUS1 * MinEgySpec / pow(SphP[i].Density / a3 , GAMMA_MINUS1);
        if(SphP[i].Entropy < minent)
            SphP[i].Entropy = minent;
    }
    MPI_Allreduce(MPI_IN_PLACE, &bad, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    if(bad > 0)
        message(0, "Detected bad densities in %d particles on disc\n",bad);
}
 
void get_mean_separation(double * MeanSeparation, const double BoxSize, const int64_t * NTotalInit)
{
    int i;
    for(i = 0; i < 6; i++) {
        if(NTotalInit[i] > 0)
            MeanSeparation[i] = BoxSize / pow(NTotalInit[i], 1.0 / 3);
    }
}
 
/* Initialize the entropy variable in Pressure-Entropy Sph.
 * Initialization of the entropy variable is a little trickier in this version of SPH,
 * since we need to make sure it 'talks to' the density appropriately */
static void
setup_density_indep_entropy(const ActiveParticles * act, ForceTree * Tree, Cosmology * CP, struct sph_pred_data * sph_pred, double u_init, double a3, int BlackHoleOn, const inttime_t Ti_Current)
{
    int j;
    int stop = 0;
 
    message(0, "Converting u -> entropy, with density split sph\n");
 
    /* This gives better convergence than initializing EgyWtDensity before Density is known*/
    #pragma omp parallel for
    for(j = 0; j < SlotsManager->info[0].size; j++)
        SphP[j].EgyWtDensity = SphP[j].Density;
 
    MyFloat * olddensity = (MyFloat *)mymalloc("olddensity ", SlotsManager->info[0].size * sizeof(MyFloat));
    for(j = 0; j < 100; j++)
    {
        int i;
        /* since ICs give energies, not entropies, need to iterate get this initialized correctly */
        #pragma omp parallel for
        for(i = 0; i < SlotsManager->info[0].size; i++) {
            SphP[i].Entropy = GAMMA_MINUS1 * u_init / pow(SphP[i].EgyWtDensity / a3 , GAMMA_MINUS1);
            olddensity[i] = SphP[i].EgyWtDensity;
        }
        /* Empty kick factors as we do not move*/
        DriftKickTimes times = init_driftkicktime(Ti_Current);
        /* Update the EgyWtDensity*/
        density(act, 0, DensityIndependentSphOn(), BlackHoleOn, 0,  times, CP, sph_pred, NULL, Tree);
        if(stop)
            break;
 
        double maxdiff = 0;
        #pragma omp parallel for reduction(max: maxdiff)
        for(i = 0; i < SlotsManager->info[0].size; i++) {
            double value = fabs(SphP[i].EgyWtDensity - olddensity[i]) / SphP[i].EgyWtDensity;
            maxdiff = DMAX(maxdiff,value);
        }
        MPI_Allreduce(MPI_IN_PLACE, &maxdiff, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
 
        message(0, "iteration %d, max relative change in EgyWtDensity = %g \n", j, maxdiff);
 
        /* If maxdiff is small, do one more iteration and stop*/
        if(maxdiff < 1e-3)
            stop = 1;
 
    }
    myfree(olddensity);
}
 
void
setup_smoothinglengths(int RestartSnapNum, DomainDecomp * ddecomp, Cosmology * CP, int BlackHoleOn, double MinEgySpec, double uu_in_cgs, const inttime_t Ti_Current, const double atime, const int64_t NTotGasInit)
{
    int i;
    const double a3 = pow(atime, 3);
 
    if(RestartSnapNum >= 0)
        return;
 
    if(InitParams.InitGasTemp < 0)
        InitParams.InitGasTemp = CP->CMBTemperature / atime;
 
    const double MeanGasSeparation = PartManager->BoxSize / pow(NTotGasInit, 1.0 / 3);
 
    int64_t tot_sph, tot_bh;
    MPI_Allreduce(&SlotsManager->info[0].size, &tot_sph, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&SlotsManager->info[5].size, &tot_bh, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
 
    /* Do nothing if we are a pure DM run*/
    if(tot_sph + tot_bh == 0)
        return;
//
 
    ForceTree Tree = {0};
    /* Need moments because we use them to set Hsml*/
    force_tree_rebuild(&Tree, ddecomp, 0, 1, NULL);
 
    /* quick hack to adjust for the baryon fraction
        * only this fraction of mass is of that type.
        * this won't work for non-dm non baryon;
        * ideally each node shall have separate count of
        * ptypes of each type.
        *
        * Eventually the iteration will fix this. */
    const double massfactor = CP->OmegaBaryon / CP->Omega0;
    const double DesNumNgb = GetNumNgb(GetDensityKernelType());
 
    #pragma omp parallel for
    for(i = 0; i < PartManager->NumPart; i++)
    {
        /* These initial smoothing lengths are only used for SPH-like particles.*/
        if(P[i].Type != 0 && P[i].Type != 5)
            continue;
 
        int no = force_get_father(i, &Tree);
 
        while(10 * DesNumNgb * P[i].Mass > massfactor * Tree.Nodes[no].mom.mass)
        {
            int p = force_get_father(no, &Tree);
 
            if(p < 0)
                break;
 
            no = p;
        }
 
        P[i].Hsml =
            pow(3.0 / (4 * M_PI) * DesNumNgb * P[i].Mass / (massfactor * Tree.Nodes[no].mom.mass),
                    1.0 / 3) * Tree.Nodes[no].len;
 
        /* recover from a poor initial guess */
        if(P[i].Hsml > 500.0 * MeanGasSeparation)
            P[i].Hsml = MeanGasSeparation;
 
        if(P[i].Hsml <= 0)
            endrun(5, "Bad hsml guess: i=%d, mass = %g type %d hsml %g no %d len %d treemass %g\n",
                    i, P[i].Mass, P[i].Type, P[i].Hsml, no, Tree.Nodes[no].len, Tree.Nodes[no].mom.mass);
    }
 
    /* for clean IC with U input only, we need to iterate to find entropy */
    double u_init = (1.0 / GAMMA_MINUS1) * (BOLTZMANN / PROTONMASS) * InitParams.InitGasTemp;
    u_init /= uu_in_cgs; /* unit conversion */
 
    double molecular_weight;
    if(InitParams.InitGasTemp > 1.0e4)  /* assuming FULL ionization */
        molecular_weight = 4 / (8 - 5 * (1 - HYDROGEN_MASSFRAC));
    else                /* assuming NEUTRAL GAS */
        molecular_weight = 4 / (1 + 3 * HYDROGEN_MASSFRAC);
    u_init /= molecular_weight;
 
    if(u_init < MinEgySpec)
        u_init = MinEgySpec;
    /* snapshot already has EgyWtDensity; hope it is read in correctly.
        * (need a test on this!) */
    /*Allocate the extra SPH data for transient SPH particle properties.*/
    struct sph_pred_data sph_pred = slots_allocate_sph_pred_data(SlotsManager->info[0].size);
 
    /*At the first time step all particles should be active*/
    ActiveParticles act = {0};
    act.ActiveParticle = NULL;
    act.NumActiveParticle = PartManager->NumPart;
 
    /* Empty kick factors as we do not move*/
    DriftKickTimes times = init_driftkicktime(Ti_Current);
    density(&act, 1, 0, BlackHoleOn, 0,  times, CP, &sph_pred, NULL, &Tree);
 
    if(DensityIndependentSphOn()) {
        setup_density_indep_entropy(&act, &Tree, CP, &sph_pred, u_init, a3, Ti_Current, BlackHoleOn);
    }
    else {
        /*Initialize to initial energy*/
        #pragma omp parallel for
        for(i = 0; i < SlotsManager->info[0].size; i++)
            SphP[i].Entropy = GAMMA_MINUS1 * u_init / pow(SphP[i].Density / a3 , GAMMA_MINUS1);
    }
    slots_free_sph_pred_data(&sph_pred);
    force_tree_free(&Tree);
}