sfr_eff.c File Reference

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <omp.h>
#include "physconst.h"
#include "sfr_eff.h"
#include "cooling.h"
#include "slotsmanager.h"
#include "walltime.h"
#include "winds.h"
#include "hydra.h"
#include "forcetree.h"
#include "domain.h"

Include dependency graph for sfr_eff.c:

Go to the source code of this file.

Classes
struct	SFRParams
struct	sfr_eeqos_data

Functions
int	get_generations (void)
static struct sfr_eeqos_data	get_sfr_eeqos (struct particle_data *part, struct sph_particle_data *sph, double dtime, struct UVBG *local_uvbg, const double redshift, const double a3inv)
static void	cooling_direct (int i, const double redshift, const double a3inv, const double hubble, const struct UVBG *const GlobalUVBG)
static void	cooling_relaxed (int i, double dtime, struct UVBG *local_uvbg, const double redshift, const double a3inv, struct sfr_eeqos_data sfr_data, const struct UVBG *const GlobalUVBG)
static int	make_particle_star (int child, int parent, int placement, double Time)
static int	starformation (int i, double *localsfr, MyFloat *sm_out, MyFloat *GradRho, const double redshift, const double a3inv, const double hubble, const double GravInternal, const struct UVBG *const GlobalUVBG)
static int	quicklyastarformation (int i, const double a3inv)
static double	get_sfr_factor_due_to_selfgravity (int i, const double atime, const double a3inv, const double hubble, const double GravInternal)
static double	get_sfr_factor_due_to_h2 (int i, MyFloat *GradRho, const double atime)
static double	get_starformation_rate_full (int i, MyFloat *GradRho, struct sfr_eeqos_data sfr_data, const double atime, const double a3inv, const double hubble, const double GravInternal)
static double	get_egyeff (double redshift, double dens, struct UVBG *uvbg)
static double	find_star_mass (int i, const double avg_baryon_mass)
static int *	sfr_reserve_slots (ActiveParticles *act, int *NewStars, int NumNewStar, ForceTree *tt)
void	set_sfr_params (ParameterSet *ps)
void	cooling_and_starformation (ActiveParticles *act, double Time, double dloga, ForceTree *tree, const Cosmology *CP, MyFloat *GradRho, FILE *FdSfr)
int	sfreff_on_eeqos (const struct sph_particle_data *sph, const double a3inv)
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)
double	get_MinEgySpec (void)
void	init_cooling_and_star_formation (int CoolingOn, int StarformationOn, Cosmology *CP, const double avg_baryon_mass, const double BoxSize, const struct UnitSystem units)
int	sfr_need_to_compute_sph_grad_rho (void)
static double	ev_NH_from_GradRho (MyFloat gradrho_mag, double hsml, double rho, double include_h)

Variables
static struct SFRParams	sfr_params

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:

cooling_direct()

static void cooling_direct ( int  i, const double  redshift, const double  a3inv, const double  hubble, const struct UVBG *const  GlobalUVBG  )

static

Definition at line 444 of file sfr_eff.c.

{
    /*  the actual time-step */
    double dloga = get_dloga_for_bin(P[i].TimeBin, P[i].Ti_drift);
    double dtime = dloga / hubble;
 
    double ne = SPHP(i).Ne; /* electron abundance (gives ionization state and mean molecular weight) */
 
    const double enttou = pow(SPHP(i).Density * a3inv, GAMMA_MINUS1) / GAMMA_MINUS1;
 
    /* Current internal energy including adiabatic change*/
    double uold = SPHP(i).Entropy * enttou;
 
    struct UVBG uvbg = get_local_UVBG(redshift, GlobalUVBG, P[i].Pos, PartManager->CurrentParticleOffset);
    double lasttime = exp(loga_from_ti(P[i].Ti_drift - dti_from_timebin(P[i].TimeBin)));
    double lastred = 1/lasttime - 1;
    double unew;
    /* The particle reionized this timestep, bump the temperature to the HI reionization temperature.
     * We only do this for non-star-forming gas.*/
    if(sfr_params.HIReionTemp > 0 && uvbg.zreion > redshift && uvbg.zreion < lastred) {
        /* We assume it is fully ionized*/
        const double meanweight = 4 / (8 - 5 * (1 - HYDROGEN_MASSFRAC));
        unew = sfr_params.temp_to_u / meanweight * sfr_params.HIReionTemp;
    }
    else {
        /* mean molecular weight assuming ZERO ionization NEUTRAL GAS*/
        const double meanweight = 4.0 / (1 + 3 * HYDROGEN_MASSFRAC);
        const double MinEgySpec = sfr_params.temp_to_u/meanweight * sfr_params.MinGasTemp;
        unew = DoCooling(redshift, uold, SPHP(i).Density * a3inv, dtime, &uvbg, &ne, SPHP(i).Metallicity, MinEgySpec, P[i].HeIIIionized);
    }
 
    SPHP(i).Ne = ne;
    /* Update the entropy. This is done after synchronizing kicks and drifts, as per run.c.*/
    SPHP(i).Entropy = unew / enttou;
    /* Cooling gas is not forming stars*/
    SPHP(i).Sfr = 0;
}

References part_manager_type::CurrentParticleOffset, dloga, DoCooling(), dti_from_timebin(), GAMMA_MINUS1, get_dloga_for_bin(), get_local_UVBG(), SFRParams::HIReionTemp, HYDROGEN_MASSFRAC, loga_from_ti(), SFRParams::MinGasTemp, P, PartManager, sfr_params, SPHP, SFRParams::temp_to_u, and UVBG::zreion.

Referenced by cooling_and_starformation().

Here is the call graph for this function:

Here is the caller graph for this function:

cooling_relaxed()

static void cooling_relaxed ( int  i, double  dtime, struct UVBG *  local_uvbg, const double  redshift, const double  a3inv, struct sfr_eeqos_data  sfr_data, const struct UVBG *const  GlobalUVBG  )

static

Definition at line 606 of file sfr_eff.c.

{
    const double egyeff = sfr_params.EgySpecCold * sfr_data.cloudfrac + (1 - sfr_data.cloudfrac) * sfr_data.egyhot;
    const double Density = SPHP(i).Density;
    const double densityfac = pow(Density * a3inv, GAMMA_MINUS1) / GAMMA_MINUS1;
    double egycurrent = SPHP(i).Entropy * densityfac;
    double trelax = sfr_data.trelax;
    if(sfr_params.BHFeedbackUseTcool == 1 && P[i].BHHeated)
    {
        if(egycurrent > egyeff)
        {
            double ne = SPHP(i).Ne;
            /* In practice tcool << trelax*/
            double tcool = GetCoolingTime(redshift, egycurrent, SPHP(i).Density * a3inv, local_uvbg, &ne, SPHP(i).Metallicity);
 
            /* The point of the star-forming equation of state is to pressurize the gas. However,
             * when the gas has been heated above the equation of state it is pressurized and does not cool successfully.
             * This code uses the cooling time rather than the relaxation time.
             * This reduces the effect of black hole feedback marginally (a 5% reduction in star formation)
             * and dates from the earliest versions of this code available.
             * The main impact is on the high end of the black hole mass function: turning this off
             * removes most massive black holes. */
            if(tcool < trelax && tcool > 0)
                trelax = tcool;
        }
        P[i].BHHeated = 0;
    }
 
    SPHP(i).Entropy =  (egyeff + (egycurrent - egyeff) * exp(-dtime / trelax)) /densityfac;
}

References SFRParams::BHFeedbackUseTcool, sfr_eeqos_data::cloudfrac, sph_particle_data::Density, sfr_eeqos_data::egyhot, SFRParams::EgySpecCold, GAMMA_MINUS1, GetCoolingTime(), sph_particle_data::Metallicity, P, sfr_params, SPHP, and sfr_eeqos_data::trelax.

Referenced by get_sfr_eeqos().

Here is the call graph for this function:

Here is the caller graph for this function:

ev_NH_from_GradRho()

static double ev_NH_from_GradRho ( MyFloat  gradrho_mag, double  hsml, double  rho, double  include_h  )

static

Definition at line 959 of file sfr_eff.c.

{
    /* column density from GradRho, copied from gadget-p; what is it
     * calculating? */
    if(rho<=0)
        return 0;
    double ev_NH = 0;
    if(gradrho_mag > 0)
        ev_NH = rho*rho/gradrho_mag;
    if(include_h > 0)
        ev_NH += rho*hsml;
    return ev_NH; // *(Z/Zsolar) add metallicity dependence
}

Referenced by get_sfr_factor_due_to_h2().

Here is the caller graph for this function:

find_star_mass()

static double find_star_mass ( int  i, const double  avg_baryon_mass  )

static

Definition at line 921 of file sfr_eff.c.

{
    /*Quick Lyman Alpha always turns all of a particle into stars*/
    if(sfr_params.QuickLymanAlphaProbability > 0)
        return P[i].Mass;
 
    double mass_of_star =  avg_baryon_mass / sfr_params.Generations;
    if(mass_of_star > P[i].Mass) {
        /* if some mass has been stolen by BH, e.g */
        mass_of_star = P[i].Mass;
    }
    /* Conditions to turn the gas into a star. .
     * The mass check makes sure we never get a gas particle which is lighter
     * than the smallest star particle.
     * The Generations check (which can happen because of mass return)
     * ensures we never instantaneously enrich stars above solar. */
    if(P[i].Mass < 2 * mass_of_star  || P[i].Generation > sfr_params.Generations) {
        mass_of_star = P[i].Mass;
    }
    return mass_of_star;
}

References SFRParams::Generations, P, SFRParams::QuickLymanAlphaProbability, and sfr_params.

Referenced by get_sfr_eeqos().

Here is the caller graph for this function:

get_egyeff()

static double get_egyeff ( double  redshift, double  dens, struct UVBG *  uvbg  )

static

Definition at line 782 of file sfr_eff.c.

{
    double tsfr = sqrt(sfr_params.PhysDensThresh / (dens)) * sfr_params.MaxSfrTimescale;
    double factorEVP = pow(dens / sfr_params.PhysDensThresh, -0.8) * sfr_params.FactorEVP;
    double egyhot = sfr_params.EgySpecSN / (1 + factorEVP) + sfr_params.EgySpecCold;
 
    double ne = 0.5;
    double tcool = GetCoolingTime(redshift, egyhot, dens, uvbg, &ne, 0.0);
 
    double y = tsfr / tcool * egyhot / (sfr_params.FactorSN * sfr_params.EgySpecSN - (1 - sfr_params.FactorSN) * sfr_params.EgySpecCold);
    double x = 1 + 1 / (2 * y) - sqrt(1 / y + 1 / (4 * y * y));
    return egyhot * (1 - x) + sfr_params.EgySpecCold * x;
}

References sfr_eeqos_data::egyhot, SFRParams::EgySpecCold, SFRParams::EgySpecSN, SFRParams::FactorEVP, SFRParams::FactorSN, GetCoolingTime(), SFRParams::MaxSfrTimescale, sfr_eeqos_data::ne, SFRParams::PhysDensThresh, sfr_params, and sfr_eeqos_data::tsfr.

Referenced by init_cooling_and_star_formation(), and sfreff_on_eeqos().

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:

get_sfr_eeqos()

struct sfr_eeqos_data get_sfr_eeqos ( struct particle_data *  part, struct sph_particle_data *  sph, double  dtime, struct UVBG *  local_uvbg, const double  redshift, const double  a3inv  )

static

Definition at line 666 of file sfr_eff.c.

{
    struct sfr_eeqos_data data;
    /* Initialise data to something, just in case.*/
    data.trelax = sfr_params.MaxSfrTimescale;
    data.tsfr = sfr_params.MaxSfrTimescale;
    data.egyhot = sfr_params.EgySpecCold;
    data.cloudfrac = 0;
    data.ne = 0;
 
    /* This shall never happen, but just in case*/
    if(!sfreff_on_eeqos(sph, a3inv))
        return data;
 
    data.ne = sph->Ne;
    data.tsfr = sqrt(sfr_params.PhysDensThresh / (sph->Density * a3inv)) * sfr_params.MaxSfrTimescale;
    /*
     * gadget-p doesn't have this cap.
     * without the cap sm can be bigger than cloudmass.
    */
    if(data.tsfr < dtime && dtime > 0)
        data.tsfr = dtime;
 
    double factorEVP = pow(sph->Density * a3inv / sfr_params.PhysDensThresh, -0.8) * sfr_params.FactorEVP;
 
    data.egyhot = sfr_params.EgySpecSN / (1 + factorEVP) + sfr_params.EgySpecCold;
    data.egycold = sfr_params.EgySpecCold;
 
    double tcool = GetCoolingTime(redshift, data.egyhot, sph->Density * a3inv, local_uvbg, &data.ne, sph->Metallicity);
    double y = data.tsfr / tcool * data.egyhot / (sfr_params.FactorSN * sfr_params.EgySpecSN - (1 - sfr_params.FactorSN) * sfr_params.EgySpecCold);
 
    data.cloudfrac = 1 + 1 / (2 * y) - sqrt(1 / y + 1 / (4 * y * y));
 
    data.trelax = data.tsfr * (1 - data.cloudfrac) / data.cloudfrac / (sfr_params.FactorSN * (1 + factorEVP));
    return data;
}

References SFRParams::avg_baryon_mass, cooling_relaxed(), part_manager_type::CurrentParticleOffset, dloga, find_star_mass(), SFRParams::Generations, get_dloga_for_bin(), get_local_UVBG(), get_random_number(), get_starformation_rate_full(), METAL_YIELD, sfr_eeqos_data::ne, P, PartManager, sfr_params, slots_split_particle(), SPHP, and SFRParams::UnitSfr_in_solar_per_year.

Referenced by get_helium_neutral_fraction_sfreff(), and get_neutral_fraction_sfreff().

Here is the call graph for this function:

Here is the caller graph for this function:

get_sfr_factor_due_to_h2()

static double get_sfr_factor_due_to_h2 ( int  i, MyFloat *  GradRho, const double  atime  )

static

Definition at line 973 of file sfr_eff.c.

                                                                                     {
    /*  Krumholz & Gnedin fitting function for f_H2 as a function of local
     *  properties, from gadget-p; we return the enhancement on SFR in this
     *  function */
    double tau_fmol;
    const double a2 = atime * atime;
    double zoverzsun = SPHP(i).Metallicity/METAL_YIELD;
    double gradrho_mag = sqrt(GradRho[3*P[i].PI]*GradRho[3*P[i].PI]+GradRho[3*P[i].PI+1]*GradRho[3*P[i].PI+1]+GradRho[3*P[i].PI+2]*GradRho[3*P[i].PI+2]);
    //message(4, "GradRho %g rho %g hsml %g i %d\n", gradrho_mag, SPHP(i).Density, P[i].Hsml, i);
    tau_fmol = ev_NH_from_GradRho(gradrho_mag,P[i].Hsml,SPHP(i).Density,1) /a2;
    tau_fmol *= (0.1 + zoverzsun);
    if(tau_fmol>0) {
        tau_fmol *= 434.78*sfr_params.tau_fmol_unit;
        double y = 0.756*(1+3.1*pow(zoverzsun,0.365));
        y = log(1+0.6*y+0.01*y*y)/(0.6*tau_fmol);
        y = 1-0.75*y/(1+0.25*y);
        if(y<0) y=0;
        if(y>1) y=1;
        return y;
 
    } // if(tau_fmol>0)
    return 1.0;
}

References ev_NH_from_GradRho(), METAL_YIELD, P, sfr_params, SPHP, and SFRParams::tau_fmol_unit.

Referenced by get_starformation_rate_full().

Here is the call graph for this function:

Here is the caller graph for this function:

get_sfr_factor_due_to_selfgravity()

static double get_sfr_factor_due_to_selfgravity ( int  i, const double  atime, const double  a3inv, const double  hubble, const double  GravInternal  )

static

Definition at line 997 of file sfr_eff.c.

                                                                                                                                               {
    const double a2 = atime * atime;
    double divv = SPHP(i).DivVel / a2;
 
    divv += 3.0*hubble * a2; // hubble-flow correction
 
    if(HAS(sfr_params.StarformationCriterion, SFR_CRITERION_CONVERGENT_FLOW)) {
        if( divv>=0 ) return 0; // restrict to convergent flows (optional) //
    }
 
    double dv2abs = (divv*divv
            + (SPHP(i).CurlVel/a2)
            * (SPHP(i).CurlVel/a2)
           ); // all in physical units
    double alpha_vir = 0.2387 * dv2abs/(GravInternal * SPHP(i).Density * a3inv);
 
    double y = 1.0;
 
    if((alpha_vir < 1.0)
    || (SPHP(i).Density * a3inv > 100. * sfr_params.PhysDensThresh)
    )  {
        y = 66.7;
    } else {
        y = 0.1;
    }
    // PFH: note the latter flag is an arbitrary choice currently set
    // -by hand- to prevent runaway densities from this prescription! //
 
    if (HAS(sfr_params.StarformationCriterion, SFR_CRITERION_CONTINUOUS_CUTOFF)) {
        // continuous cutoff w alpha_vir instead of sharp (optional) //
        y *= 1.0/(1.0 + alpha_vir);
    }
    return y;
}

References HAS, SFRParams::PhysDensThresh, SFR_CRITERION_CONTINUOUS_CUTOFF, SFR_CRITERION_CONVERGENT_FLOW, sfr_params, SPHP, and SFRParams::StarformationCriterion.

Referenced by get_starformation_rate_full().

Here is the caller graph for this function:

get_starformation_rate_full()

static double get_starformation_rate_full ( int  i, MyFloat *  GradRho, struct sfr_eeqos_data  sfr_data, const double  atime, const double  a3inv, const double  hubble, const double  GravInternal  )

static

Definition at line 759 of file sfr_eff.c.

{
    if(!sfreff_on_eeqos(&SPHP(i), a3inv)) {
        return 0;
    }
 
    double cloudmass = sfr_data.cloudfrac * P[i].Mass;
 
    double rateOfSF = (1 - sfr_params.FactorSN) * cloudmass / sfr_data.tsfr;
 
    if (HAS(sfr_params.StarformationCriterion, SFR_CRITERION_MOLECULAR_H2)) {
        if(!GradRho)
            endrun(1, "GradRho not allocated but has SFR_CRITERION_MOLECULAR_H2. Should never happen!\n");
        rateOfSF *= get_sfr_factor_due_to_h2(i, GradRho, atime);
    }
    if (HAS(sfr_params.StarformationCriterion, SFR_CRITERION_SELFGRAVITY)) {
        rateOfSF *= get_sfr_factor_due_to_selfgravity(i, atime, a3inv, hubble, GravInternal);
    }
    return rateOfSF;
}

References sfr_eeqos_data::cloudfrac, endrun(), SFRParams::FactorSN, get_sfr_factor_due_to_h2(), get_sfr_factor_due_to_selfgravity(), HAS, P, SFR_CRITERION_MOLECULAR_H2, SFR_CRITERION_SELFGRAVITY, sfr_params, sfreff_on_eeqos(), SPHP, SFRParams::StarformationCriterion, and sfr_eeqos_data::tsfr.

Referenced by get_sfr_eeqos().

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:

make_particle_star()

static int make_particle_star ( int  child, int  parent, int  placement, double  Time  )

static

Definition at line 578 of file sfr_eff.c.

{
    int retflag = 2;
    if(P[parent].Type != 0)
        endrun(7772, "Only gas forms stars, what's wrong?\n");
 
    /*Store the SPH particle slot properties, as the PI may be over-written
     *in slots_convert*/
    struct sph_particle_data oldslot = SPHP(parent);
 
    /*Convert the child slot to the new type.*/
    child = slots_convert(child, 4, placement, PartManager, SlotsManager);
 
    /*Set properties*/
    STARP(child).FormationTime = Time;
    STARP(child).LastEnrichmentMyr = 0;
    STARP(child).TotalMassReturned = 0;
    STARP(child).BirthDensity = oldslot.Density;
    /*Copy metallicity*/
    STARP(child).Metallicity = oldslot.Metallicity;
    int j;
    for(j = 0; j < NMETALS; j++)
        STARP(child).Metals[j] = oldslot.Metals[j];
    return retflag;
}

References sph_particle_data::Density, endrun(), sph_particle_data::Metallicity, sph_particle_data::Metals, NMETALS, P, PartManager, slots_convert(), SlotsManager, SPHP, and STARP.

Referenced by cooling_and_starformation().

Here is the call graph for this function:

Here is the caller graph for this function:

quicklyastarformation()

static int quicklyastarformation ( int  i, const double  a3inv  )

static

Definition at line 641 of file sfr_eff.c.

{
    if(SPHP(i).Density <= sfr_params.OverDensThresh)
        return 0;
 
    const double enttou = pow(SPHP(i).Density * a3inv, GAMMA_MINUS1) / GAMMA_MINUS1;
    double unew = SPHP(i).Entropy * enttou;
 
    const double meanweight = (4 / (8 - 5 * (1 - HYDROGEN_MASSFRAC)));
    double temp = unew * meanweight / sfr_params.temp_to_u;
 
    if(temp >= sfr_params.QuickLymanAlphaTempThresh)
        return 0;
 
    if(get_random_number(P[i].ID + 1) < sfr_params.QuickLymanAlphaProbability)
        return 1;
 
    return 0;
}

References sph_particle_data::Density, GAMMA_MINUS1, get_random_number(), HYDROGEN_MASSFRAC, SFRParams::OverDensThresh, P, SFRParams::QuickLymanAlphaProbability, SFRParams::QuickLymanAlphaTempThresh, sfr_params, SPHP, and SFRParams::temp_to_u.

Referenced by cooling_and_starformation().

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:

sfr_reserve_slots()

static int * sfr_reserve_slots ( ActiveParticles *  act, int *  NewStars, int  NumNewStar, ForceTree *  tt  )

static

Definition at line 380 of file sfr_eff.c.

{
        /* SlotsManager is below Nodes and ActiveParticleList,
         * so we need to move them out of the way before we extend Nodes.
         * This is quite slow, but need not be collective and is faster than a tree rebuild.
         * Try not to do this too often.*/
        message(1, "Need %d star slots, more than %d available. Try increasing SlotsIncreaseFactor on restart.\n", SlotsManager->info[4].size, SlotsManager->info[4].maxsize);
        /*Move the NewStar array to upper memory*/
        int * new_star_tmp = NULL;
        if(NewStars) {
            new_star_tmp = (int *) mymalloc2("newstartmp", NumNewStar*sizeof(int));
            memmove(new_star_tmp, NewStars, NumNewStar * sizeof(int));
            myfree(NewStars);
        }
        /*Move the tree to upper memory*/
        struct NODE * nodes_base_tmp=NULL;
        int *Father_tmp=NULL;
        int *ActiveParticle_tmp=NULL;
        if(force_tree_allocated(tree)) {
            nodes_base_tmp = (struct NODE *) mymalloc2("nodesbasetmp", tree->numnodes * sizeof(struct NODE));
            memmove(nodes_base_tmp, tree->Nodes_base, tree->numnodes * sizeof(struct NODE));
            myfree(tree->Nodes_base);
            Father_tmp = (int *) mymalloc2("Father_tmp", PartManager->MaxPart * sizeof(int));
            memmove(Father_tmp, tree->Father, PartManager->MaxPart * sizeof(int));
            myfree(tree->Father);
        }
        if(act->ActiveParticle) {
            ActiveParticle_tmp = (int *) mymalloc2("ActiveParticle_tmp", act->NumActiveParticle * sizeof(int));
            memmove(ActiveParticle_tmp, act->ActiveParticle, act->NumActiveParticle * sizeof(int));
            myfree(act->ActiveParticle);
        }
        /*Now we can extend the slots! */
        int64_t atleast[6];
        int64_t i;
        for(i = 0; i < 6; i++)
            atleast[i] = SlotsManager->info[i].maxsize;
        atleast[4] += NumNewStar;
        slots_reserve(1, atleast, SlotsManager);
 
        /*And now we need our memory back in the right place*/
        if(ActiveParticle_tmp) {
            act->ActiveParticle = (int *) mymalloc("ActiveParticle", sizeof(int)*(act->NumActiveParticle + PartManager->MaxPart - PartManager->NumPart));
            memmove(act->ActiveParticle, ActiveParticle_tmp, act->NumActiveParticle * sizeof(int));
            myfree(ActiveParticle_tmp);
        }
        if(force_tree_allocated(tree)) {
            tree->Father = (int *) mymalloc("Father", PartManager->MaxPart * sizeof(int));
            memmove(tree->Father, Father_tmp, PartManager->MaxPart * sizeof(int));
            myfree(Father_tmp);
            tree->Nodes_base = (struct NODE *) mymalloc("Nodes_base", tree->numnodes * sizeof(struct NODE));
            memmove(tree->Nodes_base, nodes_base_tmp, tree->numnodes * sizeof(struct NODE));
            myfree(nodes_base_tmp);
            /*Don't forget to update the Node pointer as well as Node_base!*/
            tree->Nodes = tree->Nodes_base - tree->firstnode;
        }
        if(new_star_tmp) {
            NewStars = (int *) mymalloc("NewStars", NumNewStar*sizeof(int));
            memmove(NewStars, new_star_tmp, NumNewStar * sizeof(int));
            myfree(new_star_tmp);
        }
        return NewStars;
}

References ActiveParticles::ActiveParticle, ForceTree::Father, ForceTree::firstnode, force_tree_allocated(), slots_manager_type::info, part_manager_type::MaxPart, slot_info::maxsize, message(), myfree, mymalloc, mymalloc2, ForceTree::Nodes, ForceTree::Nodes_base, ActiveParticles::NumActiveParticle, ForceTree::numnodes, part_manager_type::NumPart, PartManager, slot_info::size, slots_reserve(), and SlotsManager.

Referenced by cooling_and_starformation().

Here is the call graph for this function:

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:

starformation()

static int starformation ( int  i, double *  localsfr, MyFloat *  sm_out, MyFloat *  GradRho, const double  redshift, const double  a3inv, const double  hubble, const double  GravInternal, const struct UVBG *const  GlobalUVBG  )

static

Definition at line 666 of file sfr_eff.c.

{
    /*  the proper time-step */
    double dloga = get_dloga_for_bin(P[i].TimeBin, P[i].Ti_drift);
    double dtime = dloga / hubble;
    int newstar = -1;
 
    struct UVBG uvbg = get_local_UVBG(redshift, GlobalUVBG, P[i].Pos, PartManager->CurrentParticleOffset);
 
    struct sfr_eeqos_data sfr_data = get_sfr_eeqos(&P[i], &SPHP(i), dtime, &uvbg, redshift, a3inv);
 
    double atime = 1/(1+redshift);
    double smr = get_starformation_rate_full(i, GradRho, sfr_data, atime, a3inv, hubble, GravInternal);
 
    double sm = smr * dtime;
 
    *sm_out = sm;
    double p = sm / P[i].Mass;
 
    /* convert to Solar per Year.*/
    SPHP(i).Sfr = smr * sfr_params.UnitSfr_in_solar_per_year;
    SPHP(i).Ne = sfr_data.ne;
    *localsfr += SPHP(i).Sfr;
 
    const double w = get_random_number(P[i].ID);
    const double frac = (1 - exp(-p));
    SPHP(i).Metallicity += w * METAL_YIELD * frac / sfr_params.Generations;
 
    /* upon start-up, we need to protect against dloga ==0 */
    if(dloga > 0 && P[i].TimeBin)
        cooling_relaxed(i, dtime, &uvbg, redshift, a3inv, sfr_data, GlobalUVBG);
 
    double mass_of_star = find_star_mass(i, sfr_params.avg_baryon_mass);
    double prob = P[i].Mass / mass_of_star * (1 - exp(-p));
 
    int form_star = (get_random_number(P[i].ID + 1) < prob);
    if(form_star) {
        /* ok, make a star */
        newstar = i;
        /* If we get a fraction of the mass we need to create
         * a new particle for the star and remove mass from i.*/
        if(P[i].Mass >= 1.1 * mass_of_star)
            newstar = slots_split_particle(i, mass_of_star, PartManager);
    }
 
    /* Add the rest of the metals if we didn't form a star.
     * If we did form a star, add winds to the star-forming particle
     * that formed it if it is still around*/
    if(!form_star || newstar != i) {
        SPHP(i).Metallicity += (1-w) * METAL_YIELD * frac / sfr_params.Generations;
    }
    return newstar;
}

Referenced by cooling_and_starformation().

Here is the caller graph for this function:

Variable Documentation

sfr_params

struct SFRParams sfr_params

static

Referenced by cooling_and_starformation(), cooling_direct(), cooling_relaxed(), find_star_mass(), get_egyeff(), get_generations(), get_helium_neutral_fraction_sfreff(), get_MinEgySpec(), get_neutral_fraction_sfreff(), get_sfr_eeqos(), get_sfr_factor_due_to_h2(), get_sfr_factor_due_to_selfgravity(), get_starformation_rate_full(), init_cooling_and_star_formation(), quicklyastarformation(), set_sfr_params(), sfr_need_to_compute_sph_grad_rho(), and sfreff_on_eeqos().