timestep.h File Reference

#include "utils/paramset.h"
#include "timebinmgr.h"
#include "timefac.h"

Include dependency graph for timestep.h:

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

Go to the source code of this file.

Classes
struct	ActiveParticles
struct	DriftKickTimes

Typedefs
typedef struct ActiveParticles	ActiveParticles

Functions
int	rebuild_activelist (ActiveParticles *act, const DriftKickTimes *const times, int NumCurrentTiStep, const double Time)
void	free_activelist (ActiveParticles *act)
double	get_atime (const inttime_t Ti_Current)
int	find_timesteps (const ActiveParticles *act, DriftKickTimes *times, const double atime, int FastParticleType, const Cosmology *CP, const double asmth, const int isFirstTimeStep)
void	apply_half_kick (const ActiveParticles *act, Cosmology *CP, DriftKickTimes *times, const double atime, const double MinEgySpec)
void	apply_PM_half_kick (Cosmology *CP, DriftKickTimes *times)
int	is_timebin_active (int i, inttime_t current)
inttime_t	find_next_kick (inttime_t Ti_Current, int minTimeBin)
inttime_t	init_timebins (double TimeInit)
void	update_lastactive_drift (DriftKickTimes *times)
DriftKickTimes	init_driftkicktime (inttime_t Ti_Current)
int	is_PM_timestep (const DriftKickTimes *const times)
void	set_timestep_params (ParameterSet *ps)

Typedef Documentation

ActiveParticles

typedef struct ActiveParticles ActiveParticles

Function Documentation

apply_half_kick()

void apply_half_kick	(	const ActiveParticles *	act,
		Cosmology *	CP,
		DriftKickTimes *	times,
		const double	atime,
		const double	MinEgySpec
	)

Definition at line 323 of file timestep.c.

{
    int pa, bin;
    walltime_measure("/Misc");
    double gravkick[TIMEBINS+1] = {0}, hydrokick[TIMEBINS+1] = {0};
    /* Do nothing for the first timestep when the kicks are always zero*/
    if(times->mintimebin == 0 && times->maxtimebin == 0)
        return;
    #pragma omp parallel for
    for(bin = times->mintimebin; bin <= TIMEBINS; bin++)
    {
        /* Kick the active timebins*/
        if(is_timebin_active(bin, times->Ti_Current)) {
            /* do the kick for half a step*/
            inttime_t dti = dti_from_timebin(bin);
            inttime_t newkick = times->Ti_kick[bin] + dti/2;
            /* Compute kick factors for occupied bins*/
            gravkick[bin] = get_exact_gravkick_factor(CP, times->Ti_kick[bin], newkick);
            hydrokick[bin] = get_exact_hydrokick_factor(CP, times->Ti_kick[bin], newkick);
      //      message(0, "drift %d bin %d kick: %d->%d\n", times->Ti_Current, bin, times->Ti_kick[bin], newkick);
            times->Ti_kick[bin] = newkick;
        }
    }
    /* Advance the shorter bins without particles by the minimum occupied timestep.*/
    for(bin=1; bin < times->mintimebin; bin++)
        times->Ti_kick[bin] += dti_from_timebin(times->mintimebin)/2;
    //    message(0, "drift %d bin %d kick: %d\n", times->Ti_Current, bin, times->Ti_kick[bin]);
    /* Now assign new timesteps and kick */
    #pragma omp parallel for
    for(pa = 0; pa < act->NumActiveParticle; pa++)
    {
        const int i = get_active_particle(act, pa);
        if(P[i].Swallowed || P[i].IsGarbage)
            continue;
        int bin = P[i].TimeBin;
        if(bin > TIMEBINS)
            endrun(4, "Particle %d (type %d, id %ld) had unexpected timebin %d\n", i, P[i].Type, P[i].ID, P[i].TimeBin);
#ifdef DEBUG
        if(isnan(gravkick[bin]) || gravkick[bin] == 0.)
            endrun(5, "Bad kicks %lg bin %d tik %d\n", gravkick[bin], bin, times->Ti_kick[bin]);
#endif
        inttime_t dti = dti_from_timebin(bin);
        const double dt_entr = dloga_from_dti(dti/2, times->Ti_Current);
        /*This only changes particle i, so is thread-safe.*/
        do_the_short_range_kick(i, dt_entr, gravkick[bin], hydrokick[bin], atime, MinEgySpec);
    }
    walltime_measure("/Timeline/HalfKick/Short");
}

References CP, dloga_from_dti(), do_the_short_range_kick(), dti_from_timebin(), endrun(), get_active_particle(), get_exact_gravkick_factor(), get_exact_hydrokick_factor(), is_timebin_active(), DriftKickTimes::maxtimebin, DriftKickTimes::mintimebin, ActiveParticles::NumActiveParticle, P, DriftKickTimes::Ti_Current, DriftKickTimes::Ti_kick, TIMEBINS, and walltime_measure.

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

apply_PM_half_kick()

void apply_PM_half_kick	(	Cosmology *	CP,
		DriftKickTimes *	times
	)

Definition at line 373 of file timestep.c.

{
    /*Always do a PM half-kick, because this should be called just after a PM step*/
    const inttime_t tistart = times->PM_kick;
    const inttime_t tiend =  tistart + times->PM_length / 2;
    /* Do long-range kick */
    int i;
    const double Fgravkick = get_exact_gravkick_factor(CP, tistart, tiend);
 
    #pragma omp parallel for
    for(i = 0; i < PartManager->NumPart; i++)
    {
        int j;
        if(P[i].Swallowed || P[i].IsGarbage)
            continue;
        for(j = 0; j < 3; j++)  /* do the kick */
            P[i].Vel[j] += P[i].GravPM[j] * Fgravkick;
    }
    times->PM_kick = tiend;
    walltime_measure("/Timeline/HalfKick/Long");
}

References CP, get_exact_gravkick_factor(), GravPM, part_manager_type::NumPart, P, PartManager, DriftKickTimes::PM_kick, DriftKickTimes::PM_length, and walltime_measure.

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

find_next_kick()

inttime_t find_next_kick	(	inttime_t	Ti_Current,
		int	minTimeBin
	)

This function finds the next synchronization point of the system (i.e. the earliest point of time any of the particles needs a force computation), and drifts the system to this point of time. If the system drifts over the desired time of a snapshot file, the function will drift to this moment, generate an output, and then resume the drift.

Definition at line 712 of file timestep.c.

{
    /* Current value plus the increment for the smallest active bin. */
    return Ti_Current + dti_from_timebin(minTimeBin);
}

References dti_from_timebin().

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

find_timesteps()

int find_timesteps	(	const ActiveParticles *	act,
		DriftKickTimes *	times,
		const double	atime,
		int	FastParticleType,
		const Cosmology *	CP,
		const double	asmth,
		const int	isFirstTimeStep
	)

Definition at line 176 of file timestep.c.

{
    int pa;
    inttime_t dti_min = TIMEBASE;
 
    walltime_measure("/Misc");
 
    /*Update the PM timestep size */
    const int isPM = is_PM_timestep(times);
    inttime_t dti_max = times->PM_length;
 
    if(isPM) {
        dti_max = get_PM_timestep_ti(times, atime, CP, FastParticleType, asmth);
        times->PM_length = dti_max;
        times->PM_start = times->PM_kick;
    }
 
    const double hubble = hubble_function(CP, atime);
    /* Now assign new timesteps and kick */
    if(TimestepParams.ForceEqualTimesteps) {
        int i;
        #pragma omp parallel for reduction(min:dti_min)
        for(i = 0; i < PartManager->NumPart; i++)
        {
            enum TimeStepType titype;
            /* Because we don't GC on short timesteps, there can be garbage here.
             * Avoid making it active. */
            if(P[i].IsGarbage || P[i].Swallowed)
                continue;
            inttime_t dti = get_timestep_ti(i, dti_max, times->Ti_Current, atime, hubble, &titype);
            if(dti < dti_min)
                dti_min = dti;
        }
        MPI_Allreduce(MPI_IN_PLACE, &dti_min, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
    }
 
    int64_t ntiaccel=0, nticourant=0, ntiaccrete=0, ntineighbour=0, ntihsml=0;
    int badstepsizecount = 0;
    int mTimeBin = TIMEBINS, maxTimeBin = 0;
    #pragma omp parallel for reduction(min: mTimeBin) reduction(+: badstepsizecount, ntiaccel, nticourant, ntiaccrete, ntineighbour, ntihsml) reduction(max:maxTimeBin)
    for(pa = 0; pa < act->NumActiveParticle; pa++)
    {
        const int i = get_active_particle(act, pa);
 
        if(P[i].IsGarbage || P[i].Swallowed)
            continue;
 
        enum TimeStepType titype;
        inttime_t dti;
        if(TimestepParams.ForceEqualTimesteps) {
            dti = dti_min;
        } else {
            dti = get_timestep_ti(i, dti_max, times->Ti_Current, atime, hubble, &titype);
            if(titype == TI_ACCEL)
                ntiaccel++;
            else if (titype == TI_COURANT)
                nticourant++;
            else if (titype == TI_ACCRETE)
                ntiaccrete++;
            else if (titype == TI_NEIGH)
                ntineighbour++;
            else if (titype == TI_HSML)
                ntihsml++;
        }
 
        /* make it a power 2 subdivision */
        dti = round_down_power_of_two(dti);
 
        int bin = get_timestep_bin(dti);
        if(bin < 1) {
            message(1, "Time-step of integer size %d not allowed, id = %lu, debugging info follows.\n", dti, P[i].ID);
            badstepsizecount++;
        }
        int binold = P[i].TimeBin;
 
        /* timestep wants to increase */
        if(bin > binold)
        {
            /* make sure the new step is currently active,
             * so that particles do not miss a step */
            while(!is_timebin_active(bin, times->Ti_Current) && bin > binold && bin > 1)
                bin--;
        }
        /* This moves particles between time bins:
         * active particles always remain active
         * until rebuild_activelist is called
         * (after domain, on new timestep).*/
        P[i].TimeBin = bin;
        /*Find max and min*/
        if(bin < mTimeBin)
            mTimeBin = bin;
        if(bin > maxTimeBin)
            maxTimeBin = bin;
    }
 
    MPI_Allreduce(MPI_IN_PLACE, &badstepsizecount, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(MPI_IN_PLACE, &mTimeBin, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
    MPI_Allreduce(MPI_IN_PLACE, &maxTimeBin, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
 
    MPI_Allreduce(MPI_IN_PLACE, &ntiaccel, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(MPI_IN_PLACE, &nticourant, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(MPI_IN_PLACE, &ntihsml, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(MPI_IN_PLACE, &ntiaccrete, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(MPI_IN_PLACE, &ntineighbour, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
 
    /* Ensure that the PM timestep is not longer than the longest tree timestep;
     * this prevents particles in the longest timestep being active and moving into a higher bin
     * between PM timesteps, thus skipping the PM step entirely.*/
    if(isPM && times->PM_length > dti_from_timebin(maxTimeBin))
        times->PM_length = dti_from_timebin(maxTimeBin);
    message(0, "PM timebin: %x (dloga: %g Max: %g). Criteria: Accel: %ld Soundspeed: %ld DivVel: %ld Accrete: %ld Neighbour: %ld\n",
            times->PM_length, dloga_from_dti(times->PM_length, times->Ti_Current), TimestepParams.MaxSizeTimestep,
            ntiaccel, nticourant, ntihsml, ntiaccrete, ntineighbour);
 
    /* BH particles have their timesteps set by a timestep limiter.
     * On the first timestep this is not effective because all the particles have zero timestep.
     * So on the first timestep only set all BH particles to the smallest allowable timestep*/
    if(isFirstTimeStep) {
        #pragma omp parallel for
        for(pa = 0; pa < PartManager->NumPart; pa++)
        {
            if(P[pa].Type == 5)
                P[pa].TimeBin = mTimeBin;
        }
    }
    walltime_measure("/Timeline");
    times->mintimebin = mTimeBin;
    times->maxtimebin = maxTimeBin;
    return badstepsizecount;
}

References CP, dloga_from_dti(), dti_from_timebin(), timestep_params::ForceEqualTimesteps, get_active_particle(), get_PM_timestep_ti(), get_timestep_bin(), get_timestep_ti(), hubble_function(), is_PM_timestep(), is_timebin_active(), timestep_params::MaxSizeTimestep, DriftKickTimes::maxtimebin, message(), DriftKickTimes::mintimebin, MPI_INT64, ActiveParticles::NumActiveParticle, part_manager_type::NumPart, P, PartManager, DriftKickTimes::PM_kick, DriftKickTimes::PM_length, DriftKickTimes::PM_start, round_down_power_of_two(), TI_ACCEL, TI_ACCRETE, TI_COURANT, DriftKickTimes::Ti_Current, TI_HSML, TI_NEIGH, TIMEBASE, TIMEBINS, TimestepParams, and walltime_measure.

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

free_activelist()

void free_activelist

(

ActiveParticles *

act

)

Definition at line 795 of file timestep.c.

{
    if(act->ActiveParticle) {
        myfree(act->ActiveParticle);
    }
    event_unlisten(&EventSlotsFork, timestep_eh_slots_fork, act);
}

References ActiveParticles::ActiveParticle, event_unlisten(), EventSlotsFork, myfree, and timestep_eh_slots_fork().

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

get_atime()

double get_atime

(

const inttime_t

Ti_Current

)

Definition at line 168 of file timestep.c.

                                      {
    return exp(loga_from_ti(Ti_Current));
}

References loga_from_ti().

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

init_driftkicktime()

DriftKickTimes init_driftkicktime

(

inttime_t

Ti_Current

)

Definition at line 133 of file timestep.c.

{
    DriftKickTimes times = {0};
    times.Ti_Current = Ti_Current;
    int i;
    for(i = 0; i <= TIMEBINS; i++) {
        times.Ti_kick[i] = times.Ti_Current;
        times.Ti_lastactivedrift[i] = times.Ti_Current;
    }
    /* this makes sure the first step is a PM step. */
    times.PM_length = 0;
    times.PM_kick = times.Ti_Current;
    times.PM_start = times.Ti_Current;
    return times;
}

References DriftKickTimes::PM_kick, DriftKickTimes::PM_length, DriftKickTimes::PM_start, DriftKickTimes::Ti_Current, DriftKickTimes::Ti_kick, DriftKickTimes::Ti_lastactivedrift, and TIMEBINS.

Referenced by run(), run_gravity_test(), runfof(), setup_density_indep_entropy(), and setup_smoothinglengths().

Here is the caller graph for this function:

init_timebins()

inttime_t init_timebins

(

double

TimeInit

)

Definition at line 123 of file timestep.c.

{
    inttime_t Ti_Current = ti_from_loga(log(TimeInit));
    /*Enforce Ti_Current is initially even*/
    if(Ti_Current % 2 == 1)
        Ti_Current++;
    message(0, "Initial TimeStep at TimeInit %g Ti_Current = %d \n", TimeInit, Ti_Current);
    return Ti_Current;
}

References message(), and ti_from_loga().

Referenced by init().

Here is the call graph for this function:

Here is the caller graph for this function:

is_PM_timestep()

int is_PM_timestep

(

const DriftKickTimes *const

times

)

Definition at line 160 of file timestep.c.

{
    if(times->Ti_Current > times->PM_start + times->PM_length)
        endrun(12, "Passed end of PM step! ti=%d, PM = %d + %d\n",times->Ti_Current, times->PM_start, times->PM_length);
    return times->Ti_Current == times->PM_start + times->PM_length;
}

References endrun(), DriftKickTimes::PM_length, DriftKickTimes::PM_start, and DriftKickTimes::Ti_Current.

Referenced by find_timesteps(), print_timebin_statistics(), rebuild_activelist(), and run().

Here is the call graph for this function:

Here is the caller graph for this function:

is_timebin_active()

int is_timebin_active	(	int	i,
		inttime_t	current
	)

Definition at line 149 of file timestep.c.

                                                {
    /*Bin 0 is always active and at time 0 all bins are active*/
    if(i <= 0 || current <= 0)
        return 1;
    if(current % dti_from_timebin(i) == 0)
        return 1;
    return 0;
}

References dti_from_timebin().

Referenced by apply_half_kick(), find_timesteps(), hydro_force(), print_timebin_statistics(), rebuild_activelist(), run(), timestep_eh_slots_fork(), and update_lastactive_drift().

Here is the call graph for this function:

Here is the caller graph for this function:

rebuild_activelist()

int rebuild_activelist	(	ActiveParticles *	act,
		const DriftKickTimes *const	times,
		int	NumCurrentTiStep,
		const double	Time
	)

Definition at line 721 of file timestep.c.

{
    int i;
 
    int NumThreads = omp_get_max_threads();
    /*Since we use a static schedule, only need NumPart/NumThreads elements per thread.*/
    size_t narr = PartManager->NumPart / NumThreads + NumThreads;
 
    /*We know all particles are active on a PM timestep*/
    if(is_PM_timestep(times)) {
        act->ActiveParticle = NULL;
        act->NumActiveParticle = PartManager->NumPart;
    }
    else {
        /*Need space for more particles than we have, because of star formation*/
        act->ActiveParticle = (int *) mymalloc("ActiveParticle", narr * NumThreads * sizeof(int));
        act->NumActiveParticle = 0;
    }
 
    int * TimeBinCountType = (int *) mymalloc("TimeBinCountType", 6*(TIMEBINS+1)*NumThreads * sizeof(int));
    memset(TimeBinCountType, 0, 6 * (TIMEBINS+1) * NumThreads * sizeof(int));
 
    /*We want a lockless algorithm which preserves the ordering of the particle list.*/
    size_t *NActiveThread = ta_malloc("NActiveThread", size_t, NumThreads);
    int **ActivePartSets = ta_malloc("ActivePartSets", int *, NumThreads);
    gadget_setup_thread_arrays(act->ActiveParticle, ActivePartSets, NActiveThread, narr, NumThreads);
 
    /* We enforce schedule static to imply monotonic, ensure that each thread executes on contiguous particles
     * and ensure no thread gets more than narr particles.*/
    size_t schedsz = PartManager->NumPart / NumThreads + 1;
    #pragma omp parallel for schedule(static, schedsz)
    for(i = 0; i < PartManager->NumPart; i++)
    {
        const int bin = P[i].TimeBin;
        const int tid = omp_get_thread_num();
        if(P[i].IsGarbage || P[i].Swallowed)
            continue;
        /* when we are in PM, all particles must have been synced. */
        if (P[i].Ti_drift != times->Ti_Current) {
            endrun(5, "Particle %d type %d has drift time %x not ti_current %x!",i, P[i].Type, P[i].Ti_drift, times->Ti_Current);
        }
 
        if(act->ActiveParticle && is_timebin_active(bin, times->Ti_Current))
        {
            /* Store this particle in the ActiveSet for this thread*/
            ActivePartSets[tid][NActiveThread[tid]] = i;
            NActiveThread[tid]++;
        }
        TimeBinCountType[(TIMEBINS + 1) * (6* tid + P[i].Type) + bin] ++;
    }
    if(act->ActiveParticle) {
        /*Now we want a merge step for the ActiveParticle list.*/
        act->NumActiveParticle = gadget_compact_thread_arrays(act->ActiveParticle, ActivePartSets, NActiveThread, NumThreads);
    }
    ta_free(ActivePartSets);
    ta_free(NActiveThread);
 
    /*Print statistics for this time bin*/
    print_timebin_statistics(times, NumCurrentTiStep, TimeBinCountType, Time);
    myfree(TimeBinCountType);
 
    /* Shrink the ActiveParticle array. We still need extra space for star formation,
     * but we do not need space for the known-inactive particles*/
    if(act->ActiveParticle) {
        act->ActiveParticle = (int *) myrealloc(act->ActiveParticle, sizeof(int)*(act->NumActiveParticle + PartManager->MaxPart - PartManager->NumPart));
        act->MaxActiveParticle = act->NumActiveParticle + PartManager->MaxPart - PartManager->NumPart;
        /* listen to the slots events such that we can set timebin of new particles */
    }
    event_listen(&EventSlotsFork, timestep_eh_slots_fork, act);
    walltime_measure("/Timeline/Active");
 
    return 0;
}

References ActiveParticles::ActiveParticle, endrun(), event_listen(), EventSlotsFork, gadget_compact_thread_arrays(), gadget_setup_thread_arrays(), is_PM_timestep(), is_timebin_active(), ActiveParticles::MaxActiveParticle, part_manager_type::MaxPart, myfree, mymalloc, myrealloc, ActiveParticles::NumActiveParticle, part_manager_type::NumPart, P, PartManager, print_timebin_statistics(), ta_free, ta_malloc, DriftKickTimes::Ti_Current, TIMEBINS, timestep_eh_slots_fork(), and walltime_measure.

Referenced by run(), and run_gravity_test().

Here is the call graph for this function:

Here is the caller graph for this function:

set_timestep_params()

void set_timestep_params

(

ParameterSet *

ps

)

Definition at line 51 of file timestep.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        TimestepParams.ErrTolIntAccuracy = param_get_double(ps, "ErrTolIntAccuracy");
        TimestepParams.MaxGasVel = param_get_double(ps, "MaxGasVel");
        TimestepParams.MaxSizeTimestep = param_get_double(ps, "MaxSizeTimestep");
 
        TimestepParams.MinSizeTimestep = param_get_double(ps, "MinSizeTimestep");
        TimestepParams.ForceEqualTimesteps = param_get_int(ps, "ForceEqualTimesteps");
        TimestepParams.MaxRMSDisplacementFac = param_get_double(ps, "MaxRMSDisplacementFac");
        TimestepParams.CourantFac = param_get_double(ps, "CourantFac");
    }
    MPI_Bcast(&TimestepParams, sizeof(struct timestep_params), MPI_BYTE, 0, MPI_COMM_WORLD);
}

References timestep_params::CourantFac, timestep_params::ErrTolIntAccuracy, timestep_params::ForceEqualTimesteps, timestep_params::MaxGasVel, timestep_params::MaxRMSDisplacementFac, timestep_params::MaxSizeTimestep, timestep_params::MinSizeTimestep, param_get_double(), param_get_int(), ThisTask, and TimestepParams.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

update_lastactive_drift()

void update_lastactive_drift

(

DriftKickTimes *

times

)

Definition at line 309 of file timestep.c.

{
    int bin;
    #pragma omp parallel for
    for(bin = 0; bin <= TIMEBINS; bin++)
    {
        /* Update active timestep even if no particles in it*/
        if(is_timebin_active(bin, times->Ti_Current))
            times->Ti_lastactivedrift[bin] = times->Ti_Current;
    }
}

References is_timebin_active(), DriftKickTimes::Ti_Current, DriftKickTimes::Ti_lastactivedrift, and TIMEBINS.

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function: