init.h File Reference

#include "domain.h"
#include "utils/paramset.h"
#include "timebinmgr.h"
#include "petaio.h"

Include dependency graph for init.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
inttime_t	init (int RestartSnapNum, const char *OutputDir, struct header_data *header, Cosmology *CP)
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)
void	check_density_entropy (Cosmology *CP, const double MinEgySpec, const double atime)
void	set_init_params (ParameterSet *ps)
void	init_timeline (int RestartSnapNum, double TimeMax, const struct header_data *header, const int SnapshotWithFOF)

Function Documentation

check_density_entropy()

void check_density_entropy	(	Cosmology *	CP,
		const double	MinEgySpec,
		const double	atime
	)

Definition at line 344 of file init.c.

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

References CP, GAMMA_MINUS1, GRAVITY, HUBBLE, Cosmology::HubbleParam, slots_manager_type::info, message(), Cosmology::OmegaBaryon, slot_info::size, SlotsManager, and SphP.

Referenced by begrun().

Here is the call graph for this function:

Here is the caller graph for this function:

init()

inttime_t init	(	int	RestartSnapNum,
		const char *	OutputDir,
		struct header_data *	header,
		Cosmology *	CP
	)

This function reads the initial conditions, allocates storage for the particle data, validates and initialises the particle data.

Definition at line 86 of file init.c.

{
    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;
}

References BHP, part_manager_type::BoxSize, check_omega(), check_positions(), check_smoothing_length(), CP, domain_test_id_uniqueness(), endrun(), fof_init(), get_generations(), get_mean_separation(), Cosmology::GravInternal, gravshort_set_softenings(), HYDROGEN_MASSFRAC, init_alloc_particle_slot_memory(), init_timebins(), InitParams, header_data::MassTable, MPIU_Barrier, NMETALS, header_data::NTotalInit, part_manager_type::NumPart, P, init_params::PartAllocFactor, PartManager, PEANO(), petaio_read_snapshot(), SlotsManager, SPHP, header_data::TimeSnapshot, and walltime_measure.

Referenced by begrun().

Here is the call graph for this function:

Here is the caller graph for this function:

init_timeline()

void init_timeline	(	int	RestartSnapNum,
		double	TimeMax,
		const struct header_data *	header,
		const int	SnapshotWithFOF
	)

Definition at line 56 of file init.c.

{
    /*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);
        }
    }
}

References message(), round_down_power_of_two(), setup_sync_points(), ti_from_loga(), TIMEBASE, header_data::TimeIC, and header_data::TimeSnapshot.

Referenced by begrun().

Here is the call graph for this function:

Here is the caller graph for this function:

set_init_params()

void set_init_params

(

ParameterSet *

ps

)

Definition at line 45 of file init.c.

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

References init_params::InitGasTemp, InitParams, param_get_double(), init_params::PartAllocFactor, and ThisTask.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

setup_smoothinglengths()

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
	)

This function is used to find an initial smoothing length and initial entropy for each SPH particle. Entropies are set using the initial gas temperature. It guarantees that the number of neighbours will be between desired_ngb-MAXDEV and desired_ngb+MAXDEV. For simplicity, a first guess of the smoothing length is provided to the function density(), which will then iterate if needed to find the right smoothing length.

Definition at line 439 of file init.c.

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

References ActiveParticles::ActiveParticle, BOLTZMANN, part_manager_type::BoxSize, Cosmology::CMBTemperature, CP, density(), DensityIndependentSphOn(), endrun(), force_get_father(), force_tree_free(), force_tree_rebuild(), GAMMA_MINUS1, GetDensityKernelType(), GetNumNgb(), HYDROGEN_MASSFRAC, slots_manager_type::info, init_driftkicktime(), init_params::InitGasTemp, InitParams, NODE::len, NODE::mass, NODE::mom, MPI_INT64, ForceTree::Nodes, ActiveParticles::NumActiveParticle, part_manager_type::NumPart, Cosmology::Omega0, Cosmology::OmegaBaryon, P, PartManager, PROTONMASS, setup_density_indep_entropy(), slot_info::size, slots_allocate_sph_pred_data(), slots_free_sph_pred_data(), SlotsManager, and SphP.

Referenced by begrun().

Here is the call graph for this function:

Here is the caller graph for this function: