sfr_eff.h File Reference

#include "forcetree.h"
#include "utils/paramset.h"
#include "timestep.h"
#include "partmanager.h"
#include "slotsmanager.h"
#include "physconst.h"

Include dependency graph for sfr_eff.h:

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

Go to the source code of this file.

Macros
#define	METAL_YIELD   0.02

Enumerations
enum	StarformationCriterion {   SFR_CRITERION_DENSITY = 1 , SFR_CRITERION_MOLECULAR_H2 = 3 , SFR_CRITERION_SELFGRAVITY = 5 , SFR_CRITERION_CONVERGENT_FLOW = 13 ,   SFR_CRITERION_CONTINUOUS_CUTOFF = 21 }

Functions
void	set_sfr_params (ParameterSet *ps)
void	init_cooling_and_star_formation (int CoolingOn, int StarformationOn, Cosmology *CP, const double avg_baryon_mass, const double BoxSize, const struct UnitSystem units)
void	cooling_and_starformation (ActiveParticles *act, double Time, double dloga, ForceTree *tree, const Cosmology *CP, MyFloat *GradRho, FILE *FdSfr)
double	get_neutral_fraction_sfreff (double redshift, double hubble, struct particle_data *partdata, struct sph_particle_data *sphdata)
double	get_helium_neutral_fraction_sfreff (int ion, double redshift, double hubble, struct particle_data *partdata, struct sph_particle_data *sphdata)
int	sfr_need_to_compute_sph_grad_rho (void)
int	get_generations (void)
int	sfreff_on_eeqos (const struct sph_particle_data *sph, const double a3inv)
double	get_MinEgySpec (void)

Macro Definition Documentation

METAL_YIELD

#define METAL_YIELD   0.02

effective metal yield for star formation

Definition at line 11 of file sfr_eff.h.

Enumeration Type Documentation

StarformationCriterion

enum StarformationCriterion

Enumerator
SFR_CRITERION_DENSITY
SFR_CRITERION_MOLECULAR_H2
SFR_CRITERION_SELFGRAVITY
SFR_CRITERION_CONVERGENT_FLOW
SFR_CRITERION_CONTINUOUS_CUTOFF

Definition at line 16 of file sfr_eff.h.

                            {
    SFR_CRITERION_DENSITY = 1,
    SFR_CRITERION_MOLECULAR_H2 = 3, /* 2 + 1 */
    SFR_CRITERION_SELFGRAVITY = 5,  /* 4 + 1 */
    /* below are additional flags in SELFGRAVITY */
    SFR_CRITERION_CONVERGENT_FLOW = 13, /* 8 + 4 + 1 */
    SFR_CRITERION_CONTINUOUS_CUTOFF= 21, /* 16 + 4 + 1 */
};

Function Documentation

cooling_and_starformation()

void cooling_and_starformation	(	ActiveParticles *	act,
		double	Time,
		double	dloga,
		ForceTree *	tree,
		const Cosmology *	CP,
		MyFloat *	GradRho,
		FILE *	FdSfr
	)

Definition at line 167 of file sfr_eff.c.

{
    const int nthreads = omp_get_max_threads();
    /*This is a queue for the new stars and their parents, so we can reallocate the slots after the main cooling loop.*/
    int * NewStars = NULL;
    int * NewParents = NULL;
    int64_t NumNewStar = 0, NumMaybeWind = 0;
    int * MaybeWind = NULL;
    MyFloat * StellarMass = NULL;
 
    size_t *nqthrsfr = ta_malloc("nqthrsfr", size_t, nthreads);
    int **thrqueuesfr = ta_malloc("thrqueuesfr", int *, nthreads);
    int **thrqueueparent = ta_malloc("thrqueueparent", int *, nthreads);
    size_t *nqthrwind = NULL;
    int **thrqueuewind = NULL;
 
    /*Need to capture this so that when NumActiveParticle increases during the loop
     * we don't add extra loop iterations on particles with invalid slots.*/
    const int nactive = act->NumActiveParticle;
    const double a3inv = 1./(Time * Time * Time);
    const double hubble = hubble_function(CP, Time);
 
    if(sfr_params.StarformationOn) {
        NewStars = (int *) mymalloc("NewStars", nactive * sizeof(int) * nthreads);
        gadget_setup_thread_arrays(NewStars, thrqueuesfr, nqthrsfr, nactive, nthreads);
        NewParents = (int *) mymalloc2("NewParents", nactive * sizeof(int) * nthreads);
        gadget_setup_thread_arrays(NewParents, thrqueueparent, nqthrsfr, nactive, nthreads);
    }
 
    if(sfr_params.WindOn && winds_are_subgrid()) {
        nqthrwind = ta_malloc("nqthrwind", size_t, nthreads);
        thrqueuewind = ta_malloc("thrqueuewind", int *, nthreads);
        StellarMass = (MyFloat *) mymalloc("StellarMass", SlotsManager->info[0].size * sizeof(MyFloat));
        MaybeWind = (int *) mymalloc("MaybeWind", nactive * sizeof(int) * nthreads);
        gadget_setup_thread_arrays(MaybeWind, thrqueuewind, nqthrwind, nactive, nthreads);
    }
 
 
    /* Get the global UVBG for this redshift. */
    const double redshift = 1./Time - 1;
    struct UVBG GlobalUVBG = get_global_UVBG(redshift);
    double sum_sm = 0, sum_mass_stars = 0, localsfr = 0;
 
    /* First decide which stars are cooling and which starforming. If star forming we add them to a list.
     * Note the dynamic scheduling: individual particles may have very different loop iteration lengths.
     * Cooling is much slower than sfr. I tried splitting it into a separate loop instead, but this was faster.*/
    #pragma omp parallel reduction(+:localsfr) reduction(+: sum_sm) reduction(+:sum_mass_stars)
    {
        int i;
        const int tid = omp_get_thread_num();
        #pragma omp for schedule(static)
        for(i=0; i < nactive; i++)
        {
            /*Use raw particle number if active_set is null, otherwise use active_set*/
            const int p_i = act->ActiveParticle ? act->ActiveParticle[i] : i;
            /* Skip non-gas or garbage particles */
            if(P[p_i].Type != 0 || P[p_i].IsGarbage || P[p_i].Mass <= 0)
                continue;
 
            int shall_we_star_form = 0;
            if(sfr_params.StarformationOn) {
                /*Reduce delaytime for wind particles.*/
                winds_evolve(p_i, a3inv, hubble);
                /* check whether we are star forming gas.*/
                if(sfr_params.QuickLymanAlphaProbability > 0)
                    shall_we_star_form = quicklyastarformation(p_i, a3inv);
                else
                    shall_we_star_form = sfreff_on_eeqos(&SPHP(p_i), a3inv);
            }
 
            if(shall_we_star_form) {
                int newstar = -1;
                MyFloat sm = 0;
                if(sfr_params.QuickLymanAlphaProbability > 0) {
                    /*New star is always the same particle as the parent for quicklya*/
                    newstar = p_i;
                    sum_sm += P[p_i].Mass;
                    sm = P[p_i].Mass;
                } else {
                    newstar = starformation(p_i, &localsfr, &sm, GradRho, redshift, a3inv, hubble, CP->GravInternal, &GlobalUVBG);
                    sum_sm += P[p_i].Mass * (1 - exp(-sm/P[p_i].Mass));
                }
                /*Add this particle to the stellar conversion queue if necessary.*/
                if(newstar >= 0) {
                    thrqueuesfr[tid][nqthrsfr[tid]] = newstar;
                    thrqueueparent[tid][nqthrsfr[tid]] = p_i;
                    nqthrsfr[tid]++;
                }
                /* Add this particle to the queue for consideration to spawn a wind.
                 * Only for subgrid winds. */
                if(nqthrwind && newstar < 0) {
                    thrqueuewind[tid][nqthrwind[tid]] = p_i;
                    StellarMass[P[p_i].PI] = sm;
                    nqthrwind[tid]++;
                }
            }
            else
                cooling_direct(p_i, redshift, a3inv, hubble, &GlobalUVBG);
        }
    }
 
    report_memory_usage("SFR");
 
    walltime_measure("/Cooling/Cooling");
 
    /* Do subgrid winds*/
    if(sfr_params.WindOn && winds_are_subgrid()) {
        NumMaybeWind = gadget_compact_thread_arrays(MaybeWind, thrqueuewind, nqthrwind, nthreads);
        winds_subgrid(MaybeWind, NumMaybeWind, Time, hubble, tree, StellarMass);
        myfree(MaybeWind);
        myfree(StellarMass);
        ta_free(thrqueuewind);
        ta_free(nqthrwind);
    }
 
    /*Merge step for the queue.*/
    if(NewStars) {
        NumNewStar = gadget_compact_thread_arrays(NewStars, thrqueuesfr, nqthrsfr, nthreads);
        int64_t NumNewParent = gadget_compact_thread_arrays(NewParents, thrqueueparent, nqthrsfr, nthreads);
        if(NumNewStar != NumNewParent)
            endrun(3,"%lu new stars, but %lu new parents!\n",NumNewStar, NumNewParent);
        /*Shrink star memory as we keep it for the wind model*/
        NewStars = (int *) myrealloc(NewStars, sizeof(int) * NumNewStar);
    }
 
    ta_free(thrqueueparent);
    ta_free(thrqueuesfr);
    ta_free(nqthrsfr);
 
 
    if(!sfr_params.StarformationOn)
        return;
 
    /*Get some empty slots for the stars*/
    int firststarslot = SlotsManager->info[4].size;
    /* We ran out of slots! We must be forming a lot of stars.
     * There are things in the way of extending the slot list, so we have to move them.
     * The code in sfr_reserve_slots is not elegant, but I cannot think of a better way.*/
    if(sfr_params.StarformationOn && (SlotsManager->info[4].size + NumNewStar >= SlotsManager->info[4].maxsize)) {
        if(NewParents)
            NewParents = (int *) myrealloc(NewParents, sizeof(int) * NumNewStar);
        NewStars = sfr_reserve_slots(act, NewStars, NumNewStar, tree);
    }
    SlotsManager->info[4].size += NumNewStar;
 
    int stars_converted=0, stars_spawned=0;
    int i;
 
    /*Now we turn the particles into stars*/
    #pragma omp parallel for schedule(static) reduction(+:stars_converted) reduction(+:stars_spawned) reduction(+:sum_mass_stars)
    for(i=0; i < NumNewStar; i++)
    {
        int child = NewStars[i];
        int parent = NewParents[i];
        make_particle_star(child, parent, firststarslot+i, Time);
        sum_mass_stars += P[child].Mass;
        if(child == parent)
            stars_converted++;
        else
            stars_spawned++;
    }
 
    /*Done with the parents*/
    myfree(NewParents);
 
    int64_t tot_spawned=0, tot_converted=0;
    sumup_large_ints(1, &stars_spawned, &tot_spawned);
    sumup_large_ints(1, &stars_converted, &tot_converted);
 
    double total_sum_mass_stars, total_sm, totsfrrate;
 
    MPI_Reduce(&localsfr, &totsfrrate, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
    MPI_Reduce(&sum_sm, &total_sm, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
    MPI_Reduce(&sum_mass_stars, &total_sum_mass_stars, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
    if(FdSfr)
    {
        double rate = 0;
 
        if(dloga > 0)
            rate = total_sm / (dloga / hubble);
 
        /* convert to solar masses per yr */
 
        double rate_in_msunperyear = rate * sfr_params.UnitSfr_in_solar_per_year;
 
        /* Format:
         * Time = current scale factor,
         * total_sm = expected change in stellar mass this timestep,
         * totsfrrate = current star formation rate in active particles in Msun/year,
         * rate_in_msunperyear = expected stellar mass formation rate in Msun/year from total_sm,
         * total_sum_mass_stars = actual mass of stars formed this timestep (discretized total_sm) */
        fprintf(FdSfr, "%g %g %g %g %g\n", Time, total_sm, totsfrrate, rate_in_msunperyear,
                total_sum_mass_stars);
        fflush(FdSfr);
    }
 
    MPIU_Barrier(MPI_COMM_WORLD);
    if(tot_spawned || tot_converted)
        message(0, "SFR: spawned %ld stars, converted %ld gas particles into stars\n", tot_spawned, tot_converted);
 
    walltime_measure("/Cooling/StarFormation");
 
    /* Now apply the wind model using the list of new stars.*/
    if(sfr_params.WindOn && !winds_are_subgrid())
        winds_and_feedback(NewStars, NumNewStar, Time, hubble, tree);
 
    myfree(NewStars);
}

References ActiveParticles::ActiveParticle, cooling_direct(), CP, dloga, endrun(), gadget_compact_thread_arrays(), gadget_setup_thread_arrays(), get_global_UVBG(), Cosmology::GravInternal, hubble_function(), slots_manager_type::info, make_particle_star(), slot_info::maxsize, message(), MPIU_Barrier, myfree, mymalloc, mymalloc2, myrealloc, ActiveParticles::NumActiveParticle, P, quicklyastarformation(), SFRParams::QuickLymanAlphaProbability, report_memory_usage, sfr_params, sfr_reserve_slots(), sfreff_on_eeqos(), slot_info::size, SlotsManager, SPHP, starformation(), SFRParams::StarformationOn, sumup_large_ints(), ta_free, ta_malloc, SFRParams::UnitSfr_in_solar_per_year, walltime_measure, SFRParams::WindOn, winds_and_feedback(), winds_are_subgrid(), winds_evolve(), and winds_subgrid().

Referenced by run(), and runfof().

Here is the call graph for this function:

Here is the caller graph for this function:

get_generations()

int get_generations

(

void

)

Definition at line 87 of file sfr_eff.c.

{
    return sfr_params.Generations;
}

References SFRParams::Generations, and sfr_params.

Referenced by init(), and petaio_read_header_internal().

Here is the caller graph for this function:

get_helium_neutral_fraction_sfreff()

double get_helium_neutral_fraction_sfreff	(	int	ion,
		double	redshift,
		double	hubble,
		struct particle_data *	partdata,
		struct sph_particle_data *	sphdata
	)

Definition at line 548 of file sfr_eff.c.

{
    const double a3inv = pow(1+redshift,3);
    double helium;
    struct UVBG GlobalUVBG = get_global_UVBG(redshift);
    struct UVBG uvbg = get_local_UVBG(redshift, &GlobalUVBG, partdata->Pos, PartManager->CurrentParticleOffset);
    double physdens = sphdata->Density * a3inv;
 
    if(sfr_params.QuickLymanAlphaProbability > 0 || !sfreff_on_eeqos(sphdata, a3inv)) {
        /*This gets the neutral fraction for standard gas*/
        double eomdensity = sphdata->Density;
        double InternalEnergy = sphdata->Entropy / GAMMA_MINUS1 * pow(eomdensity * a3inv, GAMMA_MINUS1);
        helium = GetHeliumIonFraction(ion, InternalEnergy, physdens, &uvbg, sphdata->Ne);
    }
    else {
        /* This gets the neutral fraction for gas on the star-forming equation of state.
         * This needs special handling because the cold clouds have a different neutral
         * fraction than the hot gas*/
        double dloga = get_dloga_for_bin(partdata->TimeBin, partdata->Ti_drift);
        double dtime = dloga / hubble;
        struct sfr_eeqos_data sfr_data = get_sfr_eeqos(partdata, sphdata, dtime, &uvbg, redshift, a3inv);
        double nh0cold = GetHeliumIonFraction(ion, sfr_params.EgySpecCold, physdens, &uvbg, sfr_data.ne);
        double nh0hot = GetHeliumIonFraction(ion, sfr_data.egyhot, physdens, &uvbg, sfr_data.ne);
        helium =  nh0cold * sfr_data.cloudfrac + (1-sfr_data.cloudfrac) * nh0hot;
    }
    return helium;
}

References sfr_eeqos_data::cloudfrac, part_manager_type::CurrentParticleOffset, sph_particle_data::Density, dloga, sfr_eeqos_data::egyhot, SFRParams::EgySpecCold, sph_particle_data::Entropy, GAMMA_MINUS1, get_dloga_for_bin(), get_global_UVBG(), get_local_UVBG(), get_sfr_eeqos(), GetHeliumIonFraction(), sfr_eeqos_data::ne, sph_particle_data::Ne, PartManager, particle_data::Pos, SFRParams::QuickLymanAlphaProbability, sfr_params, sfreff_on_eeqos(), particle_data::Ti_drift, and particle_data::TimeBin.

Referenced by GTHeliumIFraction(), GTHeliumIIFraction(), and GTHeliumIIIFraction().

Here is the call graph for this function:

Here is the caller graph for this function:

get_MinEgySpec()

double get_MinEgySpec

(

void

)

Definition at line 797 of file sfr_eff.c.

{
    /* mean molecular weight assuming ZERO ionization NEUTRAL GAS*/
    const double meanweight = 4.0 / (1 + 3 * HYDROGEN_MASSFRAC);
    /*Enforces a minimum internal energy in cooling. */
    return sfr_params.temp_to_u / meanweight * sfr_params.MinGasTemp;
}

References HYDROGEN_MASSFRAC, SFRParams::MinGasTemp, sfr_params, and SFRParams::temp_to_u.

Referenced by begrun(), run(), and runfof().

Here is the caller graph for this function:

get_neutral_fraction_sfreff()

double get_neutral_fraction_sfreff	(	double	redshift,
		double	hubble,
		struct particle_data *	partdata,
		struct sph_particle_data *	sphdata
	)

Definition at line 520 of file sfr_eff.c.

{
    double nh0;
    const double a3inv = pow(1+redshift,3);
    struct UVBG GlobalUVBG = get_global_UVBG(redshift);
    struct UVBG uvbg = get_local_UVBG(redshift, &GlobalUVBG, partdata->Pos, PartManager->CurrentParticleOffset);
    double physdens = sphdata->Density * a3inv;
 
    if(sfr_params.QuickLymanAlphaProbability > 0 || !sfreff_on_eeqos(sphdata, a3inv)) {
        /*This gets the neutral fraction for standard gas*/
        double eomdensity = sphdata->Density;
        double InternalEnergy = sphdata->Entropy / GAMMA_MINUS1 * pow(eomdensity * a3inv, GAMMA_MINUS1);
        nh0 = GetNeutralFraction(InternalEnergy, physdens, &uvbg, sphdata->Ne);
    }
    else {
        /* This gets the neutral fraction for gas on the star-forming equation of state.
         * This needs special handling because the cold clouds have a different neutral
         * fraction than the hot gas*/
        double dloga = get_dloga_for_bin(partdata->TimeBin, partdata->Ti_drift);
        double dtime = dloga / hubble;
        struct sfr_eeqos_data sfr_data = get_sfr_eeqos(partdata, sphdata, dtime, &uvbg, redshift, a3inv);
        double nh0cold = GetNeutralFraction(sfr_params.EgySpecCold, physdens, &uvbg, sfr_data.ne);
        double nh0hot = GetNeutralFraction(sfr_data.egyhot, physdens, &uvbg, sfr_data.ne);
        nh0 =  nh0cold * sfr_data.cloudfrac + (1-sfr_data.cloudfrac) * nh0hot;
    }
    return nh0;
}

References sfr_eeqos_data::cloudfrac, part_manager_type::CurrentParticleOffset, sph_particle_data::Density, dloga, sfr_eeqos_data::egyhot, SFRParams::EgySpecCold, sph_particle_data::Entropy, GAMMA_MINUS1, get_dloga_for_bin(), get_global_UVBG(), get_local_UVBG(), get_sfr_eeqos(), GetNeutralFraction(), sfr_eeqos_data::ne, sph_particle_data::Ne, PartManager, particle_data::Pos, SFRParams::QuickLymanAlphaProbability, sfr_params, sfreff_on_eeqos(), particle_data::Ti_drift, and particle_data::TimeBin.

Referenced by blackhole_accretion_ngbiter(), and GTNeutralHydrogenFraction().

Here is the call graph for this function:

Here is the caller graph for this function:

init_cooling_and_star_formation()

void init_cooling_and_star_formation	(	int	CoolingOn,
		int	StarformationOn,
		Cosmology *	CP,
		const double	avg_baryon_mass,
		const double	BoxSize,
		const struct UnitSystem	units
	)

Definition at line 805 of file sfr_eff.c.

{
    struct cooling_units coolunits;
    coolunits.CoolingOn = CoolingOn;
    coolunits.density_in_phys_cgs = units.UnitDensity_in_cgs * CP->HubbleParam * CP->HubbleParam;
    coolunits.uu_in_cgs = units.UnitInternalEnergy_in_cgs;
    coolunits.tt_in_s = units.UnitTime_in_s / CP->HubbleParam;
    /* Get mean cosmic baryon density for photoheating rate from long mean free path photons */
    coolunits.rho_crit_baryon = 3 * pow(CP->HubbleParam * HUBBLE,2) * CP->OmegaBaryon / (8 * M_PI * GRAVITY);
 
    sfr_params.temp_to_u = (1.0 / GAMMA_MINUS1) * (BOLTZMANN / PROTONMASS) / units.UnitInternalEnergy_in_cgs;
 
    sfr_params.UnitSfr_in_solar_per_year = (units.UnitMass_in_g / SOLAR_MASS) / (units.UnitTime_in_s / SEC_PER_YEAR);
 
    init_cooling(sfr_params.TreeCoolFile, sfr_params.MetalCoolFile, sfr_params.ReionHistFile, coolunits, CP);
 
    if(!CoolingOn)
        return;
 
    /*Initialize the uv fluctuation table*/
    init_uvf_table(sfr_params.UVFluctuationFile, sizeof(sfr_params.UVFluctuationFile), BoxSize, units.UnitLength_in_cm);
 
    sfr_params.StarformationOn = StarformationOn;
 
    if(!StarformationOn)
        return;
 
    sfr_params.avg_baryon_mass = avg_baryon_mass;
 
    sfr_params.tau_fmol_unit = units.UnitDensity_in_cgs*CP->HubbleParam*units.UnitLength_in_cm;
    sfr_params.OverDensThresh =
        sfr_params.CritOverDensity * CP->OmegaBaryon * 3 * CP->Hubble * CP->Hubble / (8 * M_PI * CP->GravInternal);
 
    sfr_params.PhysDensThresh = sfr_params.CritPhysDensity * PROTONMASS / HYDROGEN_MASSFRAC / units.UnitDensity_in_cgs;
 
    /* mean molecular weight assuming ZERO ionization NEUTRAL GAS*/
    double meanweight = 4.0 / (1 + 3 * HYDROGEN_MASSFRAC);
    sfr_params.EgySpecCold = (sfr_params.temp_to_u/meanweight) * sfr_params.TempClouds;
 
    /* mean molecular weight assuming FULL ionization */
    meanweight = 4 / (8 - 5 * (1 - HYDROGEN_MASSFRAC));
    sfr_params.EgySpecSN = sfr_params.temp_to_u/meanweight * sfr_params.TempSupernova;
 
    if(sfr_params.PhysDensThresh == 0)
    {
        double egyhot = sfr_params.EgySpecSN / sfr_params.FactorEVP;
 
        meanweight = 4 / (8 - 5 * (1 - HYDROGEN_MASSFRAC)); /* note: assuming FULL ionization */
 
        double u4 = sfr_params.temp_to_u/meanweight * 1.0e4;
 
        double dens = 1.0e6 * 3 * CP->Hubble * CP->Hubble / (8 * M_PI * CP->GravInternal);
 
        double ne = 1.0;
 
        struct UVBG uvbg = {0};
        /*XXX: We set the threshold without metal cooling
         * and with zero ionization at z=0.
         * It probably make sense to set the parameters with
         * a metalicity dependence.
         * */
        const double tcool = GetCoolingTime(0, egyhot, dens, &uvbg, &ne, 0.0);
 
        const double coolrate = egyhot / tcool / dens;
 
        const double x = (egyhot - u4) / (egyhot - sfr_params.EgySpecCold);
 
        sfr_params.PhysDensThresh =
            x / pow(1 - x,
                    2) * (sfr_params.FactorSN * sfr_params.EgySpecSN - (1 -
                            sfr_params.FactorSN) * sfr_params.EgySpecCold) /
                        (sfr_params.MaxSfrTimescale * coolrate);
 
        message(0, "A0= %g  \n", sfr_params.FactorEVP);
        message(0, "Computed: PhysDensThresh= %g  (int units)         %g h^2 cm^-3\n", sfr_params.PhysDensThresh,
                sfr_params.PhysDensThresh / (PROTONMASS / HYDROGEN_MASSFRAC / units.UnitDensity_in_cgs));
        message(0, "EXPECTED FRACTION OF COLD GAS AT THRESHOLD = %g\n", x);
        message(0, "tcool=%g dens=%g egyhot=%g\n", tcool, dens, egyhot);
 
        dens = sfr_params.PhysDensThresh * 10;
 
        double neff;
        do
        {
            double egyeff = get_egyeff(0, dens, &uvbg);
 
            double peff = GAMMA_MINUS1 * dens * egyeff;
 
            const double fac = 1 / (log(dens * 1.025) - log(dens));
            neff = -log(peff) * fac;
 
            dens *= 1.025;
            egyeff = get_egyeff(0, dens, &uvbg);
            peff = GAMMA_MINUS1 * dens * egyeff;
 
            neff += log(peff) * fac;
        }
        while(neff > 4.0 / 3);
 
        message(0, "Run-away sets in for dens=%g\n", dens);
        message(0, "Dynamic range for quiescent star formation= %g\n", dens / sfr_params.PhysDensThresh);
 
        const double sigma = 10.0 / CP->Hubble * 1.0e-10 / pow(1.0e-3, 2);
 
        message(0, "Isotherm sheet central density: %g   z0=%g\n",
                M_PI * CP->GravInternal * sigma * sigma / (2 * GAMMA_MINUS1) / u4,
                GAMMA_MINUS1 * u4 / (2 * M_PI * CP->GravInternal * sigma));
    }
 
    if(sfr_params.WindOn) {
        init_winds(sfr_params.FactorSN, sfr_params.EgySpecSN, sfr_params.PhysDensThresh, units.UnitTime_in_s);
    }
 
}

References SFRParams::avg_baryon_mass, BOLTZMANN, cooling_units::CoolingOn, coolunits, CP, SFRParams::CritOverDensity, SFRParams::CritPhysDensity, cooling_units::density_in_phys_cgs, SFRParams::EgySpecCold, SFRParams::EgySpecSN, SFRParams::FactorEVP, SFRParams::FactorSN, GAMMA_MINUS1, get_egyeff(), GetCoolingTime(), Cosmology::GravInternal, GRAVITY, Cosmology::Hubble, HUBBLE, Cosmology::HubbleParam, HYDROGEN_MASSFRAC, init_cooling(), init_uvf_table(), init_winds(), SFRParams::MaxSfrTimescale, message(), SFRParams::MetalCoolFile, Cosmology::OmegaBaryon, SFRParams::OverDensThresh, SFRParams::PhysDensThresh, PROTONMASS, SFRParams::ReionHistFile, cooling_units::rho_crit_baryon, SEC_PER_YEAR, sfr_params, sigma, SOLAR_MASS, SFRParams::StarformationOn, SFRParams::tau_fmol_unit, SFRParams::temp_to_u, SFRParams::TempClouds, SFRParams::TempSupernova, SFRParams::TreeCoolFile, cooling_units::tt_in_s, UnitSystem::UnitDensity_in_cgs, UnitSystem::UnitInternalEnergy_in_cgs, UnitSystem::UnitLength_in_cm, UnitSystem::UnitMass_in_g, SFRParams::UnitSfr_in_solar_per_year, UnitSystem::UnitTime_in_s, cooling_units::uu_in_cgs, SFRParams::UVFluctuationFile, and SFRParams::WindOn.

Referenced by begrun().

Here is the call graph for this function:

Here is the caller graph for this function:

set_sfr_params()

void set_sfr_params

(

ParameterSet *

ps

)

Definition at line 129 of file sfr_eff.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        /*Star formation parameters*/
        sfr_params.StarformationCriterion = (enum StarformationCriterion) param_get_enum(ps, "StarformationCriterion");
        sfr_params.CritOverDensity = param_get_double(ps, "CritOverDensity");
        sfr_params.CritPhysDensity = param_get_double(ps, "CritPhysDensity");
        sfr_params.WindOn = param_get_int(ps, "WindOn");
 
        sfr_params.FactorSN = param_get_double(ps, "FactorSN");
        sfr_params.FactorEVP = param_get_double(ps, "FactorEVP");
        sfr_params.TempSupernova = param_get_double(ps, "TempSupernova");
        sfr_params.TempClouds = param_get_double(ps, "TempClouds");
        sfr_params.MaxSfrTimescale = param_get_double(ps, "MaxSfrTimescale");
        sfr_params.Generations = param_get_int(ps, "Generations");
        sfr_params.MinGasTemp = param_get_double(ps, "MinGasTemp");
        sfr_params.BHFeedbackUseTcool = param_get_int(ps, "BHFeedbackUseTcool");
        if(sfr_params.BHFeedbackUseTcool > 2 || sfr_params.BHFeedbackUseTcool < 0)
            endrun(0, "BHFeedbackUseTcool mode %d not supported\n", sfr_params.BHFeedbackUseTcool);
        /*Lyman-alpha forest parameters*/
        sfr_params.QuickLymanAlphaProbability = param_get_double(ps, "QuickLymanAlphaProbability");
        sfr_params.QuickLymanAlphaTempThresh = param_get_double(ps, "QuickLymanAlphaTempThresh");
        sfr_params.HIReionTemp = param_get_double(ps, "HIReionTemp");
 
        /* File names*/
        param_get_string2(ps, "TreeCoolFile", sfr_params.TreeCoolFile, sizeof(sfr_params.TreeCoolFile));
        param_get_string2(ps, "UVFluctuationfile", sfr_params.UVFluctuationFile, sizeof(sfr_params.UVFluctuationFile));
        param_get_string2(ps, "MetalCoolFile", sfr_params.MetalCoolFile, sizeof(sfr_params.MetalCoolFile));
        param_get_string2(ps, "ReionHistFile", sfr_params.ReionHistFile, sizeof(sfr_params.ReionHistFile));
    }
    MPI_Bcast(&sfr_params, sizeof(struct SFRParams), MPI_BYTE, 0, MPI_COMM_WORLD);
}

References SFRParams::BHFeedbackUseTcool, SFRParams::CritOverDensity, SFRParams::CritPhysDensity, endrun(), SFRParams::FactorEVP, SFRParams::FactorSN, SFRParams::Generations, SFRParams::HIReionTemp, SFRParams::MaxSfrTimescale, SFRParams::MetalCoolFile, SFRParams::MinGasTemp, param_get_double(), param_get_enum(), param_get_int(), param_get_string2(), SFRParams::QuickLymanAlphaProbability, SFRParams::QuickLymanAlphaTempThresh, SFRParams::ReionHistFile, sfr_params, SFRParams::StarformationCriterion, SFRParams::TempClouds, SFRParams::TempSupernova, ThisTask, SFRParams::TreeCoolFile, SFRParams::UVFluctuationFile, and SFRParams::WindOn.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

sfr_need_to_compute_sph_grad_rho()

int sfr_need_to_compute_sph_grad_rho

(

void

)

Definition at line 952 of file sfr_eff.c.

{
    if (HAS(sfr_params.StarformationCriterion, SFR_CRITERION_MOLECULAR_H2)) {
        return 1;
    }
    return 0;
}

References HAS, SFR_CRITERION_MOLECULAR_H2, sfr_params, and SFRParams::StarformationCriterion.

Referenced by run(), and runfof().

Here is the caller graph for this function:

sfreff_on_eeqos()

int sfreff_on_eeqos	(	const struct sph_particle_data *	sph,
		const double	a3inv
	)

Definition at line 486 of file sfr_eff.c.

{
    int flag = 0;
    /* no sfr: normal cooling*/
    if(!sfr_params.StarformationOn) {
        return 0;
    }
 
    if(sph->Density * a3inv >= sfr_params.PhysDensThresh)
        flag = 1;
 
    if(sph->Density < sfr_params.OverDensThresh)
        flag = 0;
 
    if(sph->DelayTime > 0)
        flag = 0;   /* only normal cooling for particles in the wind */
 
    /* The model from 0904.2572 makes gas not star forming if more than 0.5 dex above
     * the effective equation of state (at z=0). This in practice means black hole heated.*/
    if(flag == 1 && sfr_params.BHFeedbackUseTcool == 2) {
        //Redshift is the argument
        double redshift = pow(a3inv, 1./3.)-1;
        struct UVBG uvbg = get_global_UVBG(redshift);
        double egyeff = get_egyeff(redshift, sph->Density, &uvbg);
        const double enttou = pow(sph->Density * a3inv, GAMMA_MINUS1) / GAMMA_MINUS1;
        double unew = sph->Entropy * enttou;
        /* 0.5 dex = 10^0.5 = 3.2 */
        if(unew >= egyeff * 3.2)
            flag = 0;
    }
    return flag;
}

References SFRParams::BHFeedbackUseTcool, sph_particle_data::DelayTime, sph_particle_data::Density, sph_particle_data::Entropy, GAMMA_MINUS1, get_egyeff(), get_global_UVBG(), SFRParams::OverDensThresh, SFRParams::PhysDensThresh, sfr_params, and SFRParams::StarformationOn.

Referenced by blackhole_feedback_ngbiter(), cooling_and_starformation(), get_helium_neutral_fraction_sfreff(), get_neutral_fraction_sfreff(), and get_starformation_rate_full().

Here is the call graph for this function:

Here is the caller graph for this function: