run.h File Reference

#include "types.h"
#include "petaio.h"

Include dependency graph for run.h:

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

Go to the source code of this file.

Functions
inttime_t	begrun (const int RestartSnapNum, struct header_data *header)
void	run (const int RestartSnapNum, const inttime_t Ti_init, const struct header_data *header)
void	runtests (const int RestartSnapNum, const inttime_t Ti_init, const struct header_data *header)
void	runfof (const int RestartSnapNum, const inttime_t Ti_init, const struct header_data *header)
void	runpower (const struct header_data *header)
void	run_gravity_test (int RestartSnapNum, Cosmology *CP, const double Asmth, const int Nmesh, const int FastParticleType, const inttime_t Ti_Current, const char *OutputDir, const struct header_data *header)
int	find_last_snapshot (void)
void	set_all_global_params (ParameterSet *ps)

Function Documentation

begrun()

inttime_t begrun	(	const int	RestartSnapNum,
		struct header_data *	head
	)

This function performs the initial set-up of the simulation. First, the parameterfile is set, then routines for setting units, reading ICs/restart-files are called, auxialiary memory is allocated, etc.

Definition at line 194 of file run.c.

{
    petapm_module_init(omp_get_max_threads());
    petaio_init();
    walltime_init(&Clocks);
 
    *head = petaio_read_header(RestartSnapNum, All.OutputDir, &All.CP);
    /*Set Nmesh to triple the mean grid spacing of the dark matter by default.*/
    if(All.Nmesh  < 0)
        All.Nmesh = 3*pow(2, (int)(log(head->NTotal[1])/3./log(2)) );
    if(head->neutrinonk <= 0)
        head->neutrinonk = All.Nmesh;
 
    slots_init(All.SlotsIncreaseFactor * PartManager->MaxPart, SlotsManager);
    /* Enable the slots: stars and BHs are allocated if there are some,
     * or if some will form*/
    if(head->NTotalInit[0] > 0)
        slots_set_enabled(0, sizeof(struct sph_particle_data), SlotsManager);
    if(All.StarformationOn || head->NTotalInit[4] > 0)
        slots_set_enabled(4, sizeof(struct star_particle_data), SlotsManager);
    if(All.BlackHoleOn || head->NTotalInit[5] > 0)
        slots_set_enabled(5, sizeof(struct bh_particle_data), SlotsManager);
 
    const struct UnitSystem units = get_unitsystem(head->UnitLength_in_cm, head->UnitMass_in_g, head->UnitVelocity_in_cm_per_s);
    /* convert some physical input parameters to internal units */
    init_cosmology(&All.CP, head->TimeIC, units);
 
    check_units(&All.CP, units);
 
#ifdef DEBUG
    char * pidfile = fastpm_strdup_printf("%s/%s", All.OutputDir, "PIDs.txt");
 
    MPIU_write_pids(pidfile);
    myfree(pidfile);
#endif
 
    init_forcetree_params(All.FastParticleType);
 
    init_cooling_and_star_formation(All.CoolingOn, All.StarformationOn, &All.CP, head->MassTable[0], head->BoxSize, units);
 
    gravshort_fill_ntab(All.ShortRangeForceWindowType, All.Asmth);
 
    set_random_numbers(All.RandomSeed);
 
    if(All.LightconeOn)
        lightcone_init(&All.CP, head->TimeSnapshot, head->UnitLength_in_cm, All.OutputDir);
 
    init_timeline(RestartSnapNum, All.TimeMax, head, All.SnapshotWithFOF);
 
    /* Get the nk and do allocation. */
    if(All.CP.MassiveNuLinRespOn)
        init_neutrinos_lra(head->neutrinonk, head->TimeIC, All.TimeMax, All.CP.Omega0, &All.CP.ONu, All.CP.UnitTime_in_s, CM_PER_MPC);
 
    /* ... read initial model and initialise the times*/
    inttime_t ti_init = init(RestartSnapNum, All.OutputDir, head, &All.CP);
 
    if(RestartSnapNum < 0) {
        DomainDecomp ddecomp[1] = {0};
        domain_decompose_full(ddecomp); /* do initial domain decomposition (gives equal numbers of particles) so density() is safe*/
        /* On first run, generate smoothing lengths and set initial entropies based on CMB temperature*/
        setup_smoothinglengths(RestartSnapNum, ddecomp, &All.CP, All.BlackHoleOn, get_MinEgySpec(), units.UnitInternalEnergy_in_cgs, ti_init, head->TimeSnapshot, head->NTotalInit[0]);
        domain_free(ddecomp);
    }
    else
        /* When we restart, validate the SPH properties of the particles.
         * This also allows us to increase MinEgySpec on a restart if we choose.*/
        check_density_entropy(&All.CP, get_MinEgySpec(), head->TimeSnapshot);
 
    return ti_init;
}

References All, run_params::Asmth, run_params::BlackHoleOn, header_data::BoxSize, check_density_entropy(), check_units(), Clocks, CM_PER_MPC, run_params::CoolingOn, run_params::CP, domain_decompose_full(), domain_free(), run_params::FastParticleType, fastpm_strdup_printf(), get_MinEgySpec(), get_unitsystem(), gravshort_fill_ntab(), init(), init_cooling_and_star_formation(), init_cosmology(), init_forcetree_params(), init_neutrinos_lra(), init_timeline(), lightcone_init(), run_params::LightconeOn, Cosmology::MassiveNuLinRespOn, header_data::MassTable, part_manager_type::MaxPart, MPIU_write_pids(), myfree, header_data::neutrinonk, run_params::Nmesh, header_data::NTotal, header_data::NTotalInit, Cosmology::Omega0, Cosmology::ONu, run_params::OutputDir, PartManager, petaio_init(), petaio_read_header(), petapm_module_init(), run_params::RandomSeed, set_random_numbers(), setup_smoothinglengths(), run_params::ShortRangeForceWindowType, slots_init(), slots_set_enabled(), run_params::SlotsIncreaseFactor, SlotsManager, run_params::SnapshotWithFOF, run_params::StarformationOn, header_data::TimeIC, run_params::TimeMax, header_data::TimeSnapshot, UnitSystem::UnitInternalEnergy_in_cgs, header_data::UnitLength_in_cm, header_data::UnitMass_in_g, Cosmology::UnitTime_in_s, header_data::UnitVelocity_in_cm_per_s, and walltime_init().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

find_last_snapshot()

int find_last_snapshot

(

void

)

Definition at line 182 of file run.c.

{
    int RestartSnapNum = find_last_snapnum(All.OutputDir);
    message(0, "Last Snapshot number is %d.\n", RestartSnapNum);
    return RestartSnapNum;
}

References All, find_last_snapnum(), message(), and run_params::OutputDir.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

run()

void run	(	const int	RestartSnapNum,
		const inttime_t	ti_init,
		const struct header_data *	header
	)

This routine contains the main simulation loop that iterates over single timesteps. The loop terminates when the cpu-time limit is reached, when a ‘stop’ file is found in the output directory, or when the simulation ends because we arrived at TimeMax.

Definition at line 284 of file run.c.

{
    /*Number of timesteps performed this run*/
    int NumCurrentTiStep = 0;
    /*Is gas physics enabled?*/
    int GasEnabled = SlotsManager->info[0].enabled;
 
    HCIManager HCI_DEFAULT_MANAGER[1] = {0};
    hci_init(HCI_DEFAULT_MANAGER, All.OutputDir, All.TimeLimitCPU, All.AutoSnapshotTime, All.SnapshotWithFOF);
 
    const struct UnitSystem units = get_unitsystem(header->UnitLength_in_cm, header->UnitMass_in_g, header->UnitVelocity_in_cm_per_s);
 
    int SnapshotFileCount = RestartSnapNum;
 
    const double MinEgySpec = get_MinEgySpec();
 
    PetaPM pm = {0};
    gravpm_init_periodic(&pm, PartManager->BoxSize, All.Asmth, All.Nmesh, All.CP.GravInternal);
 
    DomainDecomp ddecomp[1] = {0};
 
    /* Stored scale factor of the next black hole seeding check*/
    double TimeNextSeedingCheck = header->TimeSnapshot;
 
    walltime_measure("/Misc");
    struct OutputFD fds;
    open_outputfiles(RestartSnapNum, &fds, All.OutputDir, All.BlackHoleOn, All.StarformationOn);
 
    write_cpu_log(NumCurrentTiStep, header->TimeSnapshot, fds.FdCPU, Clocks.ElapsedTime); /* produce some CPU usage info */
 
    DriftKickTimes times = init_driftkicktime(ti_init);
 
    double atime = get_atime(times.Ti_Current);
 
    while(1) /* main loop */
    {
        /* Find next synchronization point and the timebins active during this timestep.
         *
         * Note that on startup, P[i].TimeBin == 0 for all particles,
         * all bins except the zeroth are inactive and so we return 0 from this function.
         * This ensures we run the force calculation for the first timestep.
         */
        inttime_t Ti_Next = find_next_kick(times.Ti_Current, times.mintimebin);
        inttime_t Ti_Last = times.Ti_Current;
 
        times.Ti_Current = Ti_Next;
 
        /*Convert back to floating point time*/
        double newatime = get_atime(times.Ti_Current);
        if(newatime < atime)
            endrun(1, "Negative timestep: %g New Time: %g Old time %g!\n", newatime - atime, newatime, atime);
        atime = newatime;
 
        /* Compute the list of particles that cross a lightcone and write it to disc.*/
        if(All.LightconeOn)
            lightcone_compute(atime, PartManager->BoxSize, &All.CP, Ti_Last, Ti_Next);
 
        int is_PM = is_PM_timestep(&times);
 
        SyncPoint * next_sync; /* if we are out of planned sync points, terminate */
        SyncPoint * planned_sync; /* NULL; if the step is not a planned sync point. */
 
        next_sync = find_next_sync_point(times.Ti_Current);
        planned_sync = find_current_sync_point(times.Ti_Current);
 
        HCIAction action[1];
 
        hci_action_init(action); /* init to no action */
 
        int stop = 0;
 
        if(is_PM) {
            /* query HCI requests only on PM step; where kick and drifts are synced */
            stop = hci_query(HCI_DEFAULT_MANAGER, action);
 
            if(action->type == HCI_TERMINATE) {
                endrun(0, "Human triggered termination.\n");
            }
        }
 
        double rel_random_shift[3] = {0};
        if(NumCurrentTiStep > 0 && is_PM  && All.RandomParticleOffset > 0) {
            update_random_offset(PartManager, rel_random_shift, All.RandomParticleOffset);
        }
 
        int extradomain = is_timebin_active(times.mintimebin + All.MaxDomainTimeBinDepth, times.Ti_Current);
        /* drift and ddecomp decomposition */
        /* at first step this is a noop */
        if(extradomain || is_PM) {
            /* Sync positions of all particles */
            drift_all_particles(Ti_Last, times.Ti_Current, &All.CP, rel_random_shift);
            /* full decomposition rebuilds the domain, needs keys.*/
            domain_decompose_full(ddecomp);
        } else {
            /* FIXME: add a parameter for ddecomp_decompose_incremental */
            /* currently we drift all particles every step */
            /* If it is not a PM step, do a shorter version
             * of the ddecomp decomp which just exchanges particles.*/
            struct DriftData drift;
            drift.CP = &All.CP;
            drift.ti0 = Ti_Last;
            drift.ti1 = times.Ti_Current;
            domain_maintain(ddecomp, &drift);
        }
        update_lastactive_drift(&times);
 
 
        ActiveParticles Act = {0};
        rebuild_activelist(&Act, &times, NumCurrentTiStep, atime);
 
        set_random_numbers(All.RandomSeed + times.Ti_Current);
 
        /* Are the particle neutrinos gravitating this timestep?
         * If so we need to add them to the tree.*/
        int HybridNuTracer = hybrid_nu_tracer(&All.CP, atime);
 
        /* Collective: total number of active particles must be small enough*/
        int pairwisestep = use_pairwise_gravity(&Act, PartManager);
 
        MyFloat * GradRho = NULL;
        if(sfr_need_to_compute_sph_grad_rho())
            GradRho = (MyFloat *) mymalloc2("SPH_GradRho", sizeof(MyFloat) * 3 * SlotsManager->info[0].size);
 
        /* Need to rebuild the force tree because all TopLeaves are out of date.*/
        ForceTree Tree = {0};
        force_tree_rebuild(&Tree, ddecomp, HybridNuTracer, !pairwisestep && All.TreeGravOn, All.OutputDir);
 
        /* density() happens before gravity because it also initializes the predicted variables.
        * This ensures that prediction consistently uses the grav and hydro accel from the
        * timestep before this one, which matches Gadget-2/3. It was tested to make a small difference,
        * since prediction is only really used for artificial viscosity.
        *
        * Doing it first also means the density is up to date for
        * adaptive gravitational softenings. */
        if(GasEnabled)
        {
            /*Allocate the memory for predicted SPH data.*/
            struct sph_pred_data sph_predicted = slots_allocate_sph_pred_data(SlotsManager->info[0].size);
 
            if(All.DensityOn)
                density(&Act, 1, DensityIndependentSphOn(), All.BlackHoleOn, MinEgySpec, times, &All.CP, &sph_predicted, GradRho, &Tree);  /* computes density, and pressure */
 
            /***** update smoothing lengths in tree *****/
            force_update_hmax(Act.ActiveParticle, Act.NumActiveParticle, &Tree, ddecomp);
            /***** hydro forces *****/
            MPIU_Barrier(MPI_COMM_WORLD);
 
            /* adds hydrodynamical accelerations  and computes du/dt  */
            if(All.HydroOn)
                hydro_force(&Act, atime, &sph_predicted, MinEgySpec, times, &All.CP, &Tree);
 
            /* Scratch data cannot be used checkpoint because FOF does an exchange.*/
            slots_free_sph_pred_data(&sph_predicted);
        }
 
        /* The opening criterion for the gravtree
        * uses the *total* gravitational acceleration
        * from the last timestep, GravPM+GravAccel.
        * So we must compute GravAccel for this timestep
        * before gravpm_force() writes the PM acc. for
        * this timestep to GravPM. Note initially both
        * are zero and so the tree is opened maximally
        * on the first timestep.*/
        const double rho0 = All.CP.Omega0 * 3 * All.CP.Hubble * All.CP.Hubble / (8 * M_PI * All.CP.GravInternal);
 
        if(All.TreeGravOn) {
            /* Do a short range pairwise only step if desired*/
            if(pairwisestep) {
                struct gravshort_tree_params gtp = get_gravshort_treepar();
                grav_short_pair(&Act, &pm, &Tree, gtp.Rcut, rho0, HybridNuTracer, All.FastParticleType);
            }
            else
                grav_short_tree(&Act, &pm, &Tree, rho0, HybridNuTracer, All.FastParticleType);
        }
 
        /* We use the total gravitational acc.
        * to open the tree and total acc for the timestep.
        * Note that any of (GravAccel, GravPM,
        * HydroAccel) may change much faster than
        * the total acc.
        * We do the same as Gadget-2, but one could
        * instead use short-range tree acc. only
        * for opening angle or short-range timesteps,
        * or include hydro in the opening angle.*/
        if(is_PM)
        {
            gravpm_force(&pm, &Tree, &All.CP, atime, units.UnitLength_in_cm, All.OutputDir, header->TimeIC, All.FastParticleType);
 
            /* compute and output energy statistics if desired. */
            if(fds.FdEnergy)
                energy_statistics(fds.FdEnergy, atime, PartManager);
        }
 
        MPIU_Barrier(MPI_COMM_WORLD);
        message(0, "Forces computed.\n");
 
        /* Update velocity to Ti_Current; this synchronizes TiKick and TiDrift for the active particles
         * and sets Ti_Kick in the times structure.*/
        if(is_PM) {
            apply_PM_half_kick(&All.CP, &times);
        }
 
        /* Need a scale factor for entropy and velocity limiters*/
        apply_half_kick(&Act, &All.CP, &times, atime, MinEgySpec);
 
        /* Cooling and extra physics show up as a source term in the evolution equations.
         * Formally you can write the structure of the partial differential equations:
           dU/dt +  div(F) = S
         * where the cooling, BH and SFR are the source term S.
         * The extra physics is done after the kick, using a Strang split operator.
         * Gadget3/Arepo tries to follow the general operator splitting ansatz (often called Strang splitting).
         * Here you alternate the evolution under the operator generating the time evolution of the homogenous system (ie, without S)
         * with the operator generating the time evolution under the source function alone.
         * This means to advance the full system by dt, you first evolve dU/dt = S
         * by dt, and then dU/dt +  div(F) = 0 by dt.
         * [Actually, for second-order convergence in time, you should rather evolve S for dt/2, then the homogenous part for dt, and then S again for dt/2.]
 
         * The operator-split approach offers a number of practical advantages when the source function is stiff.
         * You can, for example, solve dU/dt = S in a robust and stable fashion with an implict solver, whereas the Gadget2 approach is severely challenged
         * and either requires an artficial clipping of the maximum allowed cooling rate, or a severe reduction of the timestep, otherwise the
         * predicted entropy due to cooling someting during the timestep can become severely wrong. Also, the source term approach can be easily
         * used to treat effectively instantaneous injections of energy (like from BHs), which is again hard to properly incorporate in the
         * time-integration approach where you want to have a "full" dU/dt all times. (Volker Springel 2020).
         */
        if(GasEnabled)
        {
            /* Do this before sfr and bh so the gas hsml always contains DesNumNgb neighbours.*/
            if(All.MetalReturnOn) {
                double AvgGasMass = All.CP.OmegaBaryon * 3 * All.CP.Hubble * All.CP.Hubble / (8 * M_PI * All.CP.GravInternal) * pow(PartManager->BoxSize, 3) / header->NTotalInit[0];
                metal_return(&Act, ddecomp, &All.CP, atime, AvgGasMass);
            }
 
            /* this will find new black hole seed halos.
             * Note: the FOF code does not know about garbage particles,
             * so ensure we do not have garbage present when we call this.
             * Also a good idea to only run it on a PM step.
             * This does not break the tree because the new black holes do not move or change mass, just type.
             * It does not matter that the velocities are half a step off because they are not used in the FoF code.*/
            if (is_PM && ((All.BlackHoleOn && atime >= TimeNextSeedingCheck) ||
                (during_helium_reionization(1/atime - 1) && need_change_helium_ionization_fraction(atime)))) {
 
                /* Seeding: builds its own tree.*/
                FOFGroups fof = fof_fof(ddecomp, 0, MPI_COMM_WORLD);
                if(All.BlackHoleOn && atime >= TimeNextSeedingCheck) {
                    fof_seed(&fof, &Act, atime, MPI_COMM_WORLD);
                    TimeNextSeedingCheck = atime * All.TimeBetweenSeedingSearch;
                }
 
                if(during_helium_reionization(1/atime - 1)) {
                    /* Helium reionization by switching on quasar bubbles*/
                    do_heiii_reionization(atime, &fof, ddecomp, &All.CP, units.UnitInternalEnergy_in_cgs, fds.FdHelium);
                }
                fof_finish(&fof);
            }
 
            if(is_PM) {
                /*Rebuild the force tree we freed in gravpm to save memory. Means might be two trees during FOF.*/
                force_tree_rebuild(&Tree, ddecomp, HybridNuTracer, 0, All.OutputDir);
            }
 
            /* Black hole accretion and feedback */
            if(All.BlackHoleOn) {
                blackhole(&Act, atime, &All.CP, &Tree, units, fds.FdBlackHoles, fds.FdBlackholeDetails);
            }
 
            /**** radiative cooling and star formation *****/
            if(All.CoolingOn)
                cooling_and_starformation(&Act, atime, get_dloga_for_bin(times.mintimebin, times.Ti_Current), &Tree, &All.CP, GradRho, fds.FdSfr);
 
        }
        /* We don't need this timestep's tree anymore.*/
        force_tree_free(&Tree);
 
        if(GradRho) {
            myfree(GradRho);
            GradRho = NULL;
        }
 
        /* If a snapshot is requested, write it.         *
         * We only attempt to output on sync points. This is the only chance where all variables are
         * synchronized in a consistent state in a K(KDDK)^mK scheme.
         */
 
        int WriteSnapshot = 0;
        int WriteFOF = 0;
 
        if(planned_sync) {
            WriteSnapshot |= planned_sync->write_snapshot;
            WriteFOF |= planned_sync->write_fof;
        }
 
        if(is_PM) { /* the if here is unnecessary but to signify checkpointing occurs only at PM steps. */
            WriteSnapshot |= action->write_snapshot;
            WriteFOF |= action->write_fof;
        }
        if(WriteSnapshot || WriteFOF) {
            /* Get a new snapshot*/
            SnapshotFileCount++;
            /* The accel may have created garbage -- collect them before writing a snapshot.
             * If we do collect, reset active list size.*/
            int compact[6] = {0};
            if(slots_gc(compact, PartManager, SlotsManager))
                Act.NumActiveParticle = PartManager->NumPart;
        }
        FOFGroups fof = {0};
        if(WriteFOF) {
            /* Compute FOF and assign GrNr so it can be written in checkpoint.*/
            fof = fof_fof(ddecomp, 1, MPI_COMM_WORLD);
        }
 
        /* WriteFOF just reminds the checkpoint code to save GroupID*/
        if(WriteSnapshot)
            write_checkpoint(SnapshotFileCount, WriteFOF, All.MetalReturnOn, atime, &All.CP, All.OutputDir, All.OutputDebugFields);
 
        /* Save FOF tables after checkpoint so that if there is a FOF save bug we have particle tables available to debug it*/
        if(WriteFOF) {
            fof_save_groups(&fof, All.OutputDir, All.FOFFileBase, SnapshotFileCount, &All.CP, atime, header->MassTable, All.MetalReturnOn, All.BlackHoleOn, MPI_COMM_WORLD);
            fof_finish(&fof);
        }
 
        write_cpu_log(NumCurrentTiStep, atime, fds.FdCPU, Clocks.ElapsedTime);    /* produce some CPU usage info */
 
        report_memory_usage("RUN");
 
        if(!next_sync || stop) {
            /* out of sync points, or a requested stop, the run has finally finished! Yay.*/
            if(action->type == HCI_TIMEOUT)
                message(0, "Stopping: not enough time for another PM step before TimeLimitCPU is reached.\n");
            break;
        }
 
        /* more steps to go. */
 
        /* assign new timesteps to the active particles,
         * now that we know they have synched TiKick and TiDrift,
         * and advance the PM timestep.*/
        const double asmth = All.Asmth * PartManager->BoxSize / All.Nmesh;
        int badtimestep = find_timesteps(&Act, &times, atime, All.FastParticleType, &All.CP, asmth, NumCurrentTiStep == 0);
        if(badtimestep) {
            message(0, "bad timestep spotted: terminating and saving snapshot.\n");
            dump_snapshot("TIMESTEP-DUMP", atime, &All.CP, All.OutputDir);
            endrun(0, "Ending due to bad timestep");
        }
 
        /* Update velocity and ti_kick to the new step, with the newly computed step size */
        apply_half_kick(&Act, &All.CP, &times, atime, MinEgySpec);
 
        if(is_PM) {
            apply_PM_half_kick(&All.CP, &times);
        }
 
        /* We can now free the active list: the new step have new active particles*/
        free_activelist(&Act);
 
        NumCurrentTiStep++;
    }
 
    close_outputfiles(&fds);
}

References ActiveParticles::ActiveParticle, All, apply_half_kick(), apply_PM_half_kick(), run_params::Asmth, run_params::AutoSnapshotTime, blackhole(), run_params::BlackHoleOn, part_manager_type::BoxSize, Clocks, close_outputfiles(), cooling_and_starformation(), run_params::CoolingOn, DriftData::CP, run_params::CP, density(), DensityIndependentSphOn(), run_params::DensityOn, do_heiii_reionization(), domain_decompose_full(), domain_maintain(), drift_all_particles(), dump_snapshot(), during_helium_reionization(), ClockTable::ElapsedTime, slot_info::enabled, endrun(), energy_statistics(), run_params::FastParticleType, OutputFD::FdBlackholeDetails, OutputFD::FdBlackHoles, OutputFD::FdCPU, OutputFD::FdEnergy, OutputFD::FdHelium, OutputFD::FdSfr, find_current_sync_point(), find_next_kick(), find_next_sync_point(), find_timesteps(), fof_finish(), fof_fof(), fof_save_groups(), fof_seed(), run_params::FOFFileBase, force_tree_free(), force_tree_rebuild(), force_update_hmax(), free_activelist(), get_atime(), get_dloga_for_bin(), get_gravshort_treepar(), get_MinEgySpec(), get_unitsystem(), grav_short_pair(), grav_short_tree(), Cosmology::GravInternal, gravpm_force(), gravpm_init_periodic(), hci_action_init(), hci_init(), hci_query(), HCI_TERMINATE, HCI_TIMEOUT, Cosmology::Hubble, hybrid_nu_tracer(), hydro_force(), run_params::HydroOn, slots_manager_type::info, init_driftkicktime(), is_PM_timestep(), is_timebin_active(), lightcone_compute(), run_params::LightconeOn, header_data::MassTable, run_params::MaxDomainTimeBinDepth, message(), metal_return(), run_params::MetalReturnOn, DriftKickTimes::mintimebin, MPIU_Barrier, myfree, mymalloc2, need_change_helium_ionization_fraction(), run_params::Nmesh, header_data::NTotalInit, ActiveParticles::NumActiveParticle, part_manager_type::NumPart, Cosmology::Omega0, Cosmology::OmegaBaryon, open_outputfiles(), run_params::OutputDebugFields, run_params::OutputDir, PartManager, run_params::RandomParticleOffset, run_params::RandomSeed, gravshort_tree_params::Rcut, rebuild_activelist(), report_memory_usage, set_random_numbers(), sfr_need_to_compute_sph_grad_rho(), slot_info::size, slots_allocate_sph_pred_data(), slots_free_sph_pred_data(), slots_gc(), SlotsManager, run_params::SnapshotWithFOF, run_params::StarformationOn, DriftData::ti0, DriftData::ti1, DriftKickTimes::Ti_Current, run_params::TimeBetweenSeedingSearch, header_data::TimeIC, run_params::TimeLimitCPU, header_data::TimeSnapshot, run_params::TreeGravOn, HCIAction::type, UnitSystem::UnitInternalEnergy_in_cgs, header_data::UnitLength_in_cm, UnitSystem::UnitLength_in_cm, header_data::UnitMass_in_g, header_data::UnitVelocity_in_cm_per_s, update_lastactive_drift(), update_random_offset(), use_pairwise_gravity(), walltime_measure, write_checkpoint(), write_cpu_log(), HCIAction::write_fof, SyncPoint::write_fof, HCIAction::write_snapshot, and SyncPoint::write_snapshot.

Referenced by main().

Here is the caller graph for this function:

run_gravity_test()

void run_gravity_test	(	int	RestartSnapNum,
		Cosmology *	CP,
		const double	Asmth,
		const int	Nmesh,
		const int	FastParticleType,
		const inttime_t	Ti_Current,
		const char *	OutputDir,
		const struct header_data *	header
	)

Definition at line 71 of file runtests.c.

{
    DomainDecomp ddecomp[1] = {0};
    domain_decompose_full(ddecomp);
 
    struct IOTable IOTable = {0};
    /* NO metals written*/
    register_io_blocks(&IOTable, 0, 0);
    register_extra_blocks(&IOTable);
 
    double (* PairAccn)[3] = (double (*) [3]) mymalloc2("PairAccns", 3*sizeof(double) * PartManager->NumPart);
 
    PetaPM pm[1] = {0};
    gravpm_init_periodic(pm, PartManager->BoxSize, Asmth, Nmesh, CP->GravInternal);
 
    int NTask;
    MPI_Comm_size(MPI_COMM_WORLD, &NTask);
    ActiveParticles Act = {0};
    DriftKickTimes times = init_driftkicktime(Ti_Current);
    rebuild_activelist(&Act, &times, 0, header->TimeSnapshot);
 
    ForceTree Tree = {0};
    force_tree_rebuild(&Tree, ddecomp, 1, 1, OutputDir);
    gravpm_force(pm, &Tree, CP, header->TimeSnapshot, header->UnitLength_in_cm, OutputDir, header->TimeIC, FastParticleType);
    force_tree_rebuild(&Tree, ddecomp, 1, 1, OutputDir);
 
    struct gravshort_tree_params origtreeacc = get_gravshort_treepar();
    struct gravshort_tree_params treeacc = origtreeacc;
    const double rho0 = CP->Omega0 * 3 * CP->Hubble * CP->Hubble / (8 * M_PI * CP->GravInternal);
    grav_short_pair(&Act, pm, &Tree, treeacc.Rcut, rho0, 0, FastParticleType);
 
    double meanacc = copy_and_mean_accn(PairAccn);
    message(0, "GravShort Pairs %s\n", GDB_format_particle(0));
    char * fname = fastpm_strdup_printf("%s/PART-pairs-%03d", OutputDir, RestartSnapNum);
 
    petaio_save_snapshot(fname, &IOTable, 0, header->TimeSnapshot, CP);
 
    treeacc.ErrTolForceAcc = 0;
    set_gravshort_treepar(treeacc);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
 
    /* This checks fully opened tree force against pair force*/
    double meanerr, maxerr;
    check_accns(&meanerr,&maxerr,PairAccn, meanacc);
    message(0, "Force error, open tree vs pairwise. max : %g mean: %g forcetol: %g\n", maxerr, meanerr, treeacc.ErrTolForceAcc);
 
    if(maxerr > 0.1)
        endrun(2, "Fully open tree force does not agree with pairwise calculation! maxerr %g > 0.1!\n", maxerr);
 
    message(0, "GravShort Tree %s\n", GDB_format_particle(0));
    fname = fastpm_strdup_printf("%s/PART-tree-open-%03d", OutputDir, RestartSnapNum);
    petaio_save_snapshot(fname, &IOTable, 0, header->TimeSnapshot, CP);
 
    /* This checks tree force against tree force with zero error (which always opens).*/
    copy_and_mean_accn(PairAccn);
 
    treeacc = origtreeacc;
    set_gravshort_treepar(treeacc);
    /* Code automatically sets the UseTreeBH parameter.*/
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
 
    fname = fastpm_strdup_printf("%s/PART-tree-%03d", OutputDir, RestartSnapNum);
    petaio_save_snapshot(fname, &IOTable, 0, header->TimeSnapshot, CP);
 
    check_accns(&meanerr,&maxerr,PairAccn, meanacc);
    message(0, "Force error, open tree vs tree.: %g mean: %g forcetol: %g\n", maxerr, meanerr, treeacc.ErrTolForceAcc);
 
    if(meanerr > treeacc.ErrTolForceAcc* 1.2)
        endrun(2, "Average force error is underestimated: %g > 1.2 * %g!\n", meanerr, treeacc.ErrTolForceAcc);
 
    copy_and_mean_accn(PairAccn);
    /* This checks the tree against a larger Rcut.*/
    treeacc.Rcut = 9.5;
    set_gravshort_treepar(treeacc);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
    fname = fastpm_strdup_printf("%s/PART-tree-rcut-%03d", OutputDir, RestartSnapNum);
    petaio_save_snapshot(fname, &IOTable, 0, header->TimeSnapshot, CP);
 
    check_accns(&meanerr,&maxerr,PairAccn, meanacc);
    message(0, "Force error, tree vs rcut.: %g mean: %g Rcut = %g\n", maxerr, meanerr, treeacc.Rcut);
 
    if(maxerr > 0.2)
        endrun(2, "Rcut decreased below desired value, error too large %g\n", maxerr);
 
    /* This checks the tree against a box with a smaller Nmesh.*/
    treeacc = origtreeacc;
    force_tree_free(&Tree);
 
    petapm_destroy(pm);
 
    gravpm_init_periodic(pm, PartManager->BoxSize, Asmth, Nmesh/2., CP->GravInternal);
    force_tree_rebuild(&Tree, ddecomp, 1, 1, OutputDir);
    gravpm_force(pm, &Tree, CP, header->TimeSnapshot, header->UnitLength_in_cm, OutputDir, header->TimeIC, FastParticleType);
    force_tree_rebuild(&Tree, ddecomp, 1, 1, OutputDir);
    set_gravshort_treepar(treeacc);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
    fname = fastpm_strdup_printf("%s/PART-tree-nmesh2-%03d", OutputDir, RestartSnapNum);
    petaio_save_snapshot(fname, &IOTable, 0, header->TimeSnapshot, CP);
 
    check_accns(&meanerr, &maxerr, PairAccn, meanacc);
    message(0, "Force error, nmesh %d vs %d: %g mean: %g \n", Nmesh, Nmesh/2, maxerr, meanerr);
 
    if(maxerr > 0.5 || meanerr > 0.05)
        endrun(2, "Nmesh sensitivity worse, something may be wrong\n");
 
    force_tree_free(&Tree);
    petapm_destroy(pm);
 
    myfree(PairAccn);
 
    destroy_io_blocks(&IOTable);
    domain_free(ddecomp);
}

References part_manager_type::BoxSize, check_accns(), copy_and_mean_accn(), CP, destroy_io_blocks(), domain_decompose_full(), domain_free(), endrun(), gravshort_tree_params::ErrTolForceAcc, fastpm_strdup_printf(), force_tree_free(), force_tree_rebuild(), GDB_format_particle(), get_gravshort_treepar(), grav_short_pair(), grav_short_tree(), Cosmology::GravInternal, gravpm_force(), gravpm_init_periodic(), Cosmology::Hubble, init_driftkicktime(), message(), myfree, mymalloc2, NTask, part_manager_type::NumPart, Cosmology::Omega0, PartManager, petaio_save_snapshot(), petapm_destroy(), gravshort_tree_params::Rcut, rebuild_activelist(), register_extra_blocks(), register_io_blocks(), set_gravshort_treepar(), header_data::TimeIC, header_data::TimeSnapshot, and header_data::UnitLength_in_cm.

Referenced by runtests().

Here is the call graph for this function:

Here is the caller graph for this function:

runfof()

void runfof	(	const int	RestartSnapNum,
		const inttime_t	Ti_init,
		const struct header_data *	header
	)

Definition at line 652 of file run.c.

{
    PetaPM pm = {0};
    gravpm_init_periodic(&pm, PartManager->BoxSize, All.Asmth, All.Nmesh, All.CP.GravInternal);
    DomainDecomp ddecomp[1] = {0};
    /* ... read in initial model */
 
    domain_decompose_full(ddecomp); /* do initial domain decomposition (gives equal numbers of particles) */
 
    DriftKickTimes times = init_driftkicktime(Ti_Current);
    /*FoF needs a tree*/
    int HybridNuGrav = hybrid_nu_tracer(&All.CP, header->TimeSnapshot);
    /* Regenerate the star formation rate for the FOF table.*/
    if(All.StarformationOn) {
        ActiveParticles Act = {0};
        Act.NumActiveParticle = PartManager->NumPart;
        MyFloat * GradRho = NULL;
        if(sfr_need_to_compute_sph_grad_rho()) {
            ForceTree gasTree = {0};
            GradRho = (MyFloat *) mymalloc2("SPH_GradRho", sizeof(MyFloat) * 3 * SlotsManager->info[0].size);
            /*Allocate the memory for predicted SPH data.*/
            struct sph_pred_data sph_predicted = slots_allocate_sph_pred_data(SlotsManager->info[0].size);
            force_tree_rebuild(&gasTree, ddecomp, HybridNuGrav, 0, All.OutputDir);
            /* computes GradRho with a treewalk. No hsml update as we are reading from a snapshot.*/
            density(&Act, 0, 0, All.BlackHoleOn, get_MinEgySpec(), times, &All.CP, &sph_predicted, GradRho, &gasTree);
            force_tree_free(&gasTree);
            slots_free_sph_pred_data(&sph_predicted);
        }
        ForceTree Tree = {0};
        cooling_and_starformation(&Act, header->TimeSnapshot, 0, &Tree, &All.CP, GradRho, NULL);
        if(GradRho)
            myfree(GradRho);
    }
    FOFGroups fof = fof_fof(ddecomp, 1, MPI_COMM_WORLD);
    fof_save_groups(&fof, All.OutputDir, All.FOFFileBase, RestartSnapNum, &All.CP, header->TimeSnapshot, header->MassTable, All.MetalReturnOn, All.BlackHoleOn, MPI_COMM_WORLD);
    fof_finish(&fof);
}

References All, run_params::Asmth, run_params::BlackHoleOn, part_manager_type::BoxSize, cooling_and_starformation(), run_params::CP, density(), domain_decompose_full(), fof_finish(), fof_fof(), fof_save_groups(), run_params::FOFFileBase, force_tree_free(), force_tree_rebuild(), get_MinEgySpec(), Cosmology::GravInternal, gravpm_init_periodic(), hybrid_nu_tracer(), slots_manager_type::info, init_driftkicktime(), header_data::MassTable, run_params::MetalReturnOn, myfree, mymalloc2, run_params::Nmesh, ActiveParticles::NumActiveParticle, part_manager_type::NumPart, run_params::OutputDir, PartManager, sfr_need_to_compute_sph_grad_rho(), slot_info::size, slots_allocate_sph_pred_data(), slots_free_sph_pred_data(), SlotsManager, run_params::StarformationOn, and header_data::TimeSnapshot.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

runpower()

void runpower

(

const struct header_data *

header

)

Definition at line 691 of file run.c.

{
    PetaPM pm = {0};
    gravpm_init_periodic(&pm, PartManager->BoxSize, All.Asmth, All.Nmesh, All.CP.GravInternal);
    DomainDecomp ddecomp[1] = {0};
    /* ... read in initial model */
    domain_decompose_full(ddecomp); /* do initial domain decomposition (gives equal numbers of particles) */
 
    /*PM needs a tree*/
    ForceTree Tree = {0};
    int HybridNuGrav = hybrid_nu_tracer(&All.CP, header->TimeSnapshot);
    force_tree_rebuild(&Tree, ddecomp, HybridNuGrav, 1, All.OutputDir);
    gravpm_force(&pm, &Tree, &All.CP, header->TimeSnapshot, header->UnitLength_in_cm, All.OutputDir, header->TimeSnapshot, All.FastParticleType);
    force_tree_free(&Tree);
}

References All, run_params::Asmth, part_manager_type::BoxSize, run_params::CP, domain_decompose_full(), run_params::FastParticleType, force_tree_free(), force_tree_rebuild(), Cosmology::GravInternal, gravpm_force(), gravpm_init_periodic(), hybrid_nu_tracer(), run_params::Nmesh, run_params::OutputDir, PartManager, header_data::TimeSnapshot, and header_data::UnitLength_in_cm.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

runtests()

void runtests	(	const int	RestartSnapNum,
		const inttime_t	Ti_init,
		const struct header_data *	header
	)

Definition at line 646 of file run.c.

{
    run_gravity_test(RestartSnapNum, &All.CP, All.Asmth, All.Nmesh, All.FastParticleType, Ti_Current, All.OutputDir, header);
}

References All, run_params::Asmth, run_params::CP, run_params::FastParticleType, run_params::Nmesh, run_params::OutputDir, and run_gravity_test().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

set_all_global_params()

void set_all_global_params

(

ParameterSet *

ps

)

Definition at line 106 of file run.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        /* Start reading the values */
        param_get_string2(ps, "OutputDir", All.OutputDir, sizeof(All.OutputDir));
        param_get_string2(ps, "FOFFileBase", All.FOFFileBase, sizeof(All.FOFFileBase));
 
        All.CP.CMBTemperature = param_get_double(ps, "CMBTemperature");
        All.CP.RadiationOn = param_get_int(ps, "RadiationOn");
        All.CP.Omega0 = param_get_double(ps, "Omega0");
        All.CP.OmegaBaryon = param_get_double(ps, "OmegaBaryon");
        All.CP.OmegaLambda = param_get_double(ps, "OmegaLambda");
        All.CP.Omega_fld = param_get_double(ps, "Omega_fld");
        if(All.CP.OmegaLambda > 0 && All.CP.Omega_fld > 0)
            endrun(0, "Cannot have OmegaLambda and Omega_fld (evolving dark energy) at the same time!\n");
        All.CP.w0_fld = param_get_double(ps,"w0_fld");
        All.CP.wa_fld = param_get_double(ps,"wa_fld");
        All.CP.Omega_ur = param_get_double(ps, "Omega_ur");
        All.CP.HubbleParam = param_get_double(ps, "HubbleParam");
 
        All.OutputDebugFields = param_get_int(ps, "OutputDebugFields");
 
        All.TimeMax = param_get_double(ps, "TimeMax");
        All.Asmth = param_get_double(ps, "Asmth");
        All.ShortRangeForceWindowType = (enum ShortRangeForceWindowType) param_get_enum(ps, "ShortRangeForceWindowType");
        All.Nmesh = param_get_int(ps, "Nmesh");
 
        All.CoolingOn = param_get_int(ps, "CoolingOn");
        All.HydroOn = param_get_int(ps, "HydroOn");
        All.DensityOn = param_get_int(ps, "DensityOn");
        All.TreeGravOn = param_get_int(ps, "TreeGravOn");
        All.LightconeOn = param_get_int(ps, "LightconeOn");
        All.FastParticleType = param_get_int(ps, "FastParticleType");
        All.PairwiseActiveFraction = param_get_double(ps, "PairwiseActiveFraction");
        All.TimeLimitCPU = param_get_double(ps, "TimeLimitCPU");
        All.AutoSnapshotTime = param_get_double(ps, "AutoSnapshotTime");
        All.TimeBetweenSeedingSearch = param_get_double(ps, "TimeBetweenSeedingSearch");
        All.RandomParticleOffset = param_get_double(ps, "RandomParticleOffset");
        /* Convert to a fraction of the box, from a fraction of a PM mesh cell*/
        All.RandomParticleOffset /= All.Nmesh;
 
        All.SlotsIncreaseFactor = param_get_double(ps, "SlotsIncreaseFactor");
 
        All.SnapshotWithFOF = param_get_int(ps, "SnapshotWithFOF");
 
        All.RandomSeed = param_get_int(ps, "RandomSeed");
 
        All.BlackHoleOn = param_get_int(ps, "BlackHoleOn");
 
        All.StarformationOn = param_get_int(ps, "StarformationOn");
        All.MetalReturnOn = param_get_int(ps, "MetalReturnOn");
        All.MaxDomainTimeBinDepth = param_get_int(ps, "MaxDomainTimeBinDepth");
 
        /*Massive neutrino parameters*/
        All.CP.MassiveNuLinRespOn = param_get_int(ps, "MassiveNuLinRespOn");
        All.CP.HybridNeutrinosOn = param_get_int(ps, "HybridNeutrinosOn");
        All.CP.MNu[0] = param_get_double(ps, "MNue");
        All.CP.MNu[1] = param_get_double(ps, "MNum");
        All.CP.MNu[2] = param_get_double(ps, "MNut");
        All.CP.HybridVcrit = param_get_double(ps, "Vcrit");
        All.CP.HybridNuPartTime = param_get_double(ps, "NuPartTime");
        if(All.CP.MassiveNuLinRespOn && !All.CP.RadiationOn)
            endrun(2, "You have enabled (kspace) massive neutrinos without radiation, but this will give an inconsistent cosmology!\n");
        /*End massive neutrino parameters*/
 
        if(All.StarformationOn != 0 && All.CoolingOn == 0)
        {
                endrun(1, "You try to use the code with star formation enabled,\n"
                          "but you did not switch on cooling.\nThis mode is not supported.\n");
        }
    }
    MPI_Bcast(&All, sizeof(All), MPI_BYTE, 0, MPI_COMM_WORLD);
}

References All, run_params::Asmth, run_params::AutoSnapshotTime, run_params::BlackHoleOn, Cosmology::CMBTemperature, run_params::CoolingOn, run_params::CP, run_params::DensityOn, endrun(), run_params::FastParticleType, run_params::FOFFileBase, Cosmology::HubbleParam, Cosmology::HybridNeutrinosOn, Cosmology::HybridNuPartTime, Cosmology::HybridVcrit, run_params::HydroOn, run_params::LightconeOn, Cosmology::MassiveNuLinRespOn, run_params::MaxDomainTimeBinDepth, run_params::MetalReturnOn, Cosmology::MNu, run_params::Nmesh, Cosmology::Omega0, Cosmology::Omega_fld, Cosmology::Omega_ur, Cosmology::OmegaBaryon, Cosmology::OmegaLambda, run_params::OutputDebugFields, run_params::OutputDir, run_params::PairwiseActiveFraction, param_get_double(), param_get_enum(), param_get_int(), param_get_string2(), Cosmology::RadiationOn, run_params::RandomParticleOffset, run_params::RandomSeed, run_params::ShortRangeForceWindowType, run_params::SlotsIncreaseFactor, run_params::SnapshotWithFOF, run_params::StarformationOn, ThisTask, run_params::TimeBetweenSeedingSearch, run_params::TimeLimitCPU, run_params::TimeMax, run_params::TreeGravOn, Cosmology::w0_fld, and Cosmology::wa_fld.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function: