winds.h File Reference

#include "forcetree.h"
#include "utils/paramset.h"

Include dependency graph for winds.h:

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

Go to the source code of this file.

Enumerations
enum	WindModel {   WIND_SUBGRID = 1 , WIND_DECOUPLE_SPH = 2 , WIND_USE_HALO = 4 , WIND_FIXED_EFFICIENCY = 8 ,   WIND_ISOTROPIC = 512 }

Functions
void	set_winds_params (ParameterSet *ps)
void	init_winds (double FactorSN, double EgySpecSN, double PhysDensThresh, double UnitTime_in_s)
void	winds_evolve (int i, double a3inv, double hubble)
void	winds_and_feedback (int *NewStars, int NumNewStars, const double Time, const double hubble, ForceTree *tree)
int	winds_make_after_sf (int i, double sm, double vdisp, double atime)
void	winds_subgrid (int *MaybeWind, int NumMaybeWind, const double Time, const double hubble, ForceTree *tree, MyFloat *StellarMasses)
int	winds_are_subgrid (void)
int	winds_is_particle_decoupled (int i)
void	winds_decoupled_hydro (int i, double atime)
int	winds_ever_decouple (void)

Enumeration Type Documentation

WindModel

enum WindModel

Enumerator
WIND_SUBGRID
WIND_DECOUPLE_SPH
WIND_USE_HALO
WIND_FIXED_EFFICIENCY
WIND_ISOTROPIC

Definition at line 12 of file winds.h.

               {
    WIND_SUBGRID = 1,/* If this is true, winds are spawned from the star forming gas.
                      * If false, they are spawned from neighbours of the star particle.*/
    WIND_DECOUPLE_SPH = 2,/* Specifies that wind particles are created temporarily decoupled from the gas dynamics */
    WIND_USE_HALO = 4,/* Wind speeds depend on the halo circular velocity*/
    WIND_FIXED_EFFICIENCY = 8, /* Winds have a fixed efficiency and thus fixed wind speed*/
    WIND_ISOTROPIC = 512, /* Has no effect: only isotropic winds are implemented*/
};

Function Documentation

init_winds()

void init_winds	(	double	FactorSN,
		double	EgySpecSN,
		double	PhysDensThresh,
		double	UnitTime_in_s
	)

Definition at line 97 of file winds.c.

{
    wind_params.WindSpeed = sqrt(2 * wind_params.WindEnergyFraction * FactorSN * EgySpecSN / (1 - FactorSN));
    /* Convert wind free travel time from Myr to internal units*/
    wind_params.MaxWindFreeTravelTime = wind_params.MaxWindFreeTravelTime * SEC_PER_MEGAYEAR / UnitTime_in_s;
    wind_params.WindFreeTravelDensThresh = wind_params.WindFreeTravelDensFac * PhysDensThresh;
    if(HAS(wind_params.WindModel, WIND_FIXED_EFFICIENCY)) {
        wind_params.WindSpeed /= sqrt(wind_params.WindEfficiency);
        message(0, "Windspeed: %g MaxDelay %g\n", wind_params.WindSpeed, wind_params.MaxWindFreeTravelTime);
    } else if(HAS(wind_params.WindModel, WIND_USE_HALO)) {
        message(0, "Reference Windspeed: %g, MaxDelay %g\n", wind_params.WindSigma0 * wind_params.WindSpeedFactor, wind_params.MaxWindFreeTravelTime);
    } else {
        /* Check for undefined wind models*/
        endrun(1, "WindModel = 0x%X is strange. This shall not happen.\n", wind_params.WindModel);
    }
}

References endrun(), HAS, WindParams::MaxWindFreeTravelTime, message(), SEC_PER_MEGAYEAR, WIND_FIXED_EFFICIENCY, wind_params, WIND_USE_HALO, WindParams::WindEfficiency, WindParams::WindEnergyFraction, WindParams::WindFreeTravelDensFac, WindParams::WindFreeTravelDensThresh, WindParams::WindModel, WindParams::WindSigma0, WindParams::WindSpeed, and WindParams::WindSpeedFactor.

Referenced by init_cooling_and_star_formation().

Here is the call graph for this function:

Here is the caller graph for this function:

set_winds_params()

void set_winds_params

(

ParameterSet *

ps

)

Definition at line 72 of file winds.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        /*Wind model parameters*/
        wind_params.WindModel = (enum WindModel) param_get_enum(ps, "WindModel");
        /* The following two are for VS08 and SH03*/
        wind_params.WindEfficiency = param_get_double(ps, "WindEfficiency");
        wind_params.WindEnergyFraction = param_get_double(ps, "WindEnergyFraction");
 
        /* The following two are for OFJT10*/
        wind_params.WindSigma0 = param_get_double(ps, "WindSigma0");
        wind_params.WindSpeedFactor = param_get_double(ps, "WindSpeedFactor");
 
        wind_params.WindThermalFactor = param_get_double(ps, "WindThermalFactor");
        wind_params.MinWindVelocity = param_get_double(ps, "MinWindVelocity");
        wind_params.MaxWindFreeTravelTime = param_get_double(ps, "MaxWindFreeTravelTime");
        wind_params.WindFreeTravelLength = param_get_double(ps, "WindFreeTravelLength");
        wind_params.WindFreeTravelDensFac = param_get_double(ps, "WindFreeTravelDensFac");
    }
    MPI_Bcast(&wind_params, sizeof(struct WindParams), MPI_BYTE, 0, MPI_COMM_WORLD);
}

References WindParams::MaxWindFreeTravelTime, WindParams::MinWindVelocity, param_get_double(), param_get_enum(), ThisTask, wind_params, WindParams::WindEfficiency, WindParams::WindEnergyFraction, WindParams::WindFreeTravelDensFac, WindParams::WindFreeTravelLength, WindParams::WindModel, WindParams::WindSigma0, WindParams::WindSpeedFactor, and WindParams::WindThermalFactor.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

winds_and_feedback()

void winds_and_feedback	(	int *	NewStars,
		int	NumNewStars,
		const double	Time,
		const double	hubble,
		ForceTree *	tree
	)

Definition at line 333 of file winds.c.

{
    /*The subgrid model does nothing here*/
    if(HAS(wind_params.WindModel, WIND_SUBGRID))
        return;
 
    if(!MPIU_Any(NumNewStars > 0, MPI_COMM_WORLD))
        return;
 
    TreeWalk tw[1] = {{0}};
    struct WindPriv priv[1];
    int i;
    winds_find_weights(tw, priv, NewStars, NumNewStars, Time, hubble, tree);
 
    for (i = 1; i < omp_get_max_threads(); i++)
        priv->nvisited[0] += priv->nvisited[i];
    priv->maxkicks = priv->nvisited[0]+2;
 
    /* Some particles may be kicked by multiple stars on the same timestep.
    * To ensure this happens only once and does not depend on the order in
    * which the loops are executed, particles are kicked by the nearest new star.
    * This struct stores all such possible kicks, and we sort it out after the treewalk.*/
    priv->kicks = (struct StarKick * ) mymalloc("StarKicks", priv->maxkicks * sizeof(struct StarKick));
    priv->nkicks = 0;
    ta_free(priv->nvisited);
 
    /* Then run feedback */
    tw->haswork = NULL;
    tw->ngbiter = (TreeWalkNgbIterFunction) sfr_wind_feedback_ngbiter;
    tw->postprocess = NULL;
    tw->reduce = NULL;
    tw->ev_label = "WIND_KICK";
    tw->Niteration = 0;
 
    treewalk_run(tw, NewStars, NumNewStars);
 
    /* Sort the possible kicks*/
    qsort_openmp(priv->kicks, priv->nkicks, sizeof(struct StarKick), cmp_by_part_id);
    /* Not parallel as the number of kicked particles should be pretty small*/
    int64_t last_part = -1;
    int64_t nkicked = 0;
    for(i = 0; i < priv->nkicks; i++) {
        /* Only do the kick for the first particle, which is the closest*/
        if(priv->kicks[i].part_index == last_part)
            continue;
        int other = priv->kicks[i].part_index;
        last_part = other;
        nkicked++;
        wind_do_kick(other, priv->kicks[i].StarKickVelocity, priv->kicks[i].StarTherm, Time);
        if(priv->kicks[i].StarKickVelocity <= 0 || !isfinite(priv->kicks[i].StarKickVelocity) || !isfinite(SPHP(other).DelayTime))
        {
            endrun(5, "Odd v: other = %d, DT = %g v = %g i = %d, nkicks %d maxkicks %d dist %g id %ld\n",
                   other, SPHP(other).DelayTime, priv->kicks[i].StarKickVelocity, i, priv->nkicks, priv->maxkicks,
                   priv->kicks[i].StarDistance, priv->kicks[i].StarID);
        }
    }
    /* Get total number of potential new stars to allocate memory.*/
    int64_t tot_newstars, tot_kicks, tot_applied;
    sumup_large_ints(1, &NumNewStars, &tot_newstars);
    MPI_Allreduce(&priv->nkicks, &tot_kicks, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&nkicked, &tot_applied, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    message(0, "Made %ld gas wind, discarded %ld kicks from %d stars\n", tot_applied, tot_kicks - tot_applied, tot_newstars);
 
    myfree(priv->kicks);
    myfree(priv->Winddata);
    walltime_measure("/Cooling/Wind");
}

References cmp_by_part_id(), endrun(), TreeWalk::ev_label, HAS, TreeWalk::haswork, WindPriv::hubble, WindPriv::kicks, WindPriv::maxkicks, message(), MPI_INT64, MPIU_Any(), myfree, mymalloc, TreeWalk::ngbiter, TreeWalk::Niteration, WindPriv::nkicks, WindPriv::nvisited, StarKick::part_index, TreeWalk::postprocess, qsort_openmp, TreeWalk::reduce, sfr_wind_feedback_ngbiter(), SPHP, StarKick::StarDistance, StarKick::StarID, StarKick::StarKickVelocity, StarKick::StarTherm, sumup_large_ints(), ta_free, WindPriv::Time, treewalk_run(), walltime_measure, wind_do_kick(), wind_params, WIND_SUBGRID, WindPriv::Winddata, WindParams::WindModel, and winds_find_weights().

Referenced by cooling_and_starformation().

Here is the call graph for this function:

Here is the caller graph for this function:

winds_are_subgrid()

int winds_are_subgrid

(

void

)

Definition at line 115 of file winds.c.

{
    if(HAS(wind_params.WindModel, WIND_SUBGRID))
        return 1;
    else
        return 0;
}

References HAS, wind_params, WIND_SUBGRID, and WindParams::WindModel.

Referenced by cooling_and_starformation().

Here is the caller graph for this function:

winds_decoupled_hydro()

void winds_decoupled_hydro	(	int	i,
		double	atime
	)

Definition at line 133 of file winds.c.

{
    int k;
    for(k = 0; k < 3; k++)
        SPHP(i).HydroAccel[k] = 0;
 
    SPHP(i).DtEntropy = 0;
 
    double windspeed = wind_params.WindSpeed * atime;
    const double fac_mu = pow(atime, 3 * (GAMMA - 1) / 2) / atime;
    windspeed *= fac_mu;
    double hsml_c = cbrt(wind_params.WindFreeTravelDensThresh /SPHP(i).Density) * atime;
    SPHP(i).MaxSignalVel = hsml_c * DMAX((2 * windspeed), SPHP(i).MaxSignalVel);
}

References DMAX, GAMMA, SPHP, wind_params, WindParams::WindFreeTravelDensThresh, and WindParams::WindSpeed.

Referenced by hydro_postprocess().

Here is the caller graph for this function:

winds_ever_decouple()

int winds_ever_decouple

(

void

)

Definition at line 192 of file winds.c.

{
    if(wind_params.MaxWindFreeTravelTime > 0)
        return 1;
    else
        return 0;
}

References WindParams::MaxWindFreeTravelTime, and wind_params.

Referenced by wind_do_kick().

Here is the caller graph for this function:

winds_evolve()

void winds_evolve	(	int	i,
		double	a3inv,
		double	hubble
	)

Definition at line 403 of file winds.c.

{
    /*Remove a wind particle from the delay mode if the (physical) density has dropped sufficiently.*/
    if(SPHP(i).DelayTime > 0 && SPHP(i).Density * a3inv < wind_params.WindFreeTravelDensThresh) {
        SPHP(i).DelayTime = 0;
    }
    /*Reduce the time until the particle can form stars again by the current timestep*/
    if(SPHP(i).DelayTime > 0) {
        /* Enforce the maximum in case of restarts*/
        if(SPHP(i).DelayTime > wind_params.MaxWindFreeTravelTime)
            SPHP(i).DelayTime = wind_params.MaxWindFreeTravelTime;
        const double dloga = get_dloga_for_bin(P[i].TimeBin, P[i].Ti_drift);
        /*  the proper time duration of the step */
        const double dtime = dloga / hubble;
        SPHP(i).DelayTime = DMAX(SPHP(i).DelayTime - dtime, 0);
    }
}

References dloga, DMAX, get_dloga_for_bin(), WindPriv::hubble, WindParams::MaxWindFreeTravelTime, P, SPHP, wind_params, and WindParams::WindFreeTravelDensThresh.

Referenced by cooling_and_starformation().

Here is the call graph for this function:

Here is the caller graph for this function:

winds_is_particle_decoupled()

int winds_is_particle_decoupled

(

int

i

)

Definition at line 124 of file winds.c.

{
    if(HAS(wind_params.WindModel, WIND_DECOUPLE_SPH)
        && P[i].Type == 0 && SPHP(i).DelayTime > 0)
            return 1;
    return 0;
}

References HAS, P, SPHP, WIND_DECOUPLE_SPH, wind_params, and WindParams::WindModel.

Referenced by add_particle_to_group(), blackhole_accretion_ngbiter(), blackhole_feedback_ngbiter(), density_ngbiter(), hydro_ngbiter(), and hydro_postprocess().

Here is the caller graph for this function:

winds_make_after_sf()

int winds_make_after_sf	(	int	i,
		double	sm,
		double	vdisp,
		double	atime
	)

Definition at line 708 of file winds.c.

{
    if(!HAS(wind_params.WindModel, WIND_SUBGRID))
        return 0;
 
    /* Get the velocity, thermal energy and efficiency of the kick*/
    double utherm = 0, vel=0, windeff = 0;
    get_wind_params(&vel, &windeff, &utherm, vdisp, atime);
 
    /* Here comes the Springel Hernquist 03 wind model */
    /* Notice that this is the mass of the gas particle after forking a star, Mass - Mass/GENERATIONS.*/
    double pw = windeff * sm / P[i].Mass;
    double prob = 1 - exp(-pw);
    if(get_random_number(P[i].ID + 2) < prob) {
        wind_do_kick(i, vel, utherm, atime);
    }
    return 0;
}

Referenced by winds_subgrid().

Here is the caller graph for this function:

winds_subgrid()

void winds_subgrid	(	int *	MaybeWind,
		int	NumMaybeWind,
		const double	Time,
		const double	hubble,
		ForceTree *	tree,
		MyFloat *	StellarMasses
	)

Definition at line 306 of file winds.c.

{
    /*The non-subgrid model does nothing here*/
    if(!HAS(wind_params.WindModel, WIND_SUBGRID))
        return;
 
    if(!MPIU_Any(NumMaybeWind > 0, MPI_COMM_WORLD))
        return;
 
    TreeWalk tw[1] = {{0}};
    struct WindPriv priv[1];
    int n;
    winds_find_weights(tw, priv, MaybeWind, NumMaybeWind, Time, hubble, tree);
    myfree(priv->nvisited);
    for(n = 0; n < NumMaybeWind; n++)
    {
        int i = MaybeWind ? MaybeWind[n] : n;
        /* Notice that StellarMasses is indexed like PI, not i!*/
        MyFloat sm = StellarMasses[P[i].PI];
        winds_make_after_sf(i, sm, WINDP(i, priv->Winddata).Vdisp, Time);
    }
    myfree(priv->Winddata);
    walltime_measure("/Cooling/Wind");
}

References HAS, WindPriv::hubble, MPIU_Any(), myfree, WindPriv::nvisited, P, WindPriv::Time, walltime_measure, wind_params, WIND_SUBGRID, WindPriv::Winddata, WindParams::WindModel, WINDP, winds_find_weights(), and winds_make_after_sf().

Referenced by cooling_and_starformation().

Here is the call graph for this function:

Here is the caller graph for this function: