cooling.h File Reference

#include "cosmology.h"

Include dependency graph for cooling.h:

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

Go to the source code of this file.

Classes
struct	UVBG
struct	cooling_units

Functions
void	init_cooling (const char *TreeCoolFile, const char *MetalCoolFile, char *reion_hist_file, struct cooling_units cu, Cosmology *CP)
void	init_uvf_table (const char *UVFluctuationFile, const int UVFlucLen, const double BoxSize, const double UnitLength_in_cm)
double	GetCoolingTime (double redshift, double u_old, double rho, struct UVBG *uvbg, double *ne_guess, double Z)
double	DoCooling (double redshift, double u_old, double rho, double dt, struct UVBG *uvbg, double *ne_guess, double Z, double MinEgySpec, int isHeIIIionized)
struct UVBG	get_global_UVBG (double redshift)
struct UVBG	get_local_UVBG (double redshift, const struct UVBG *const GlobalUVBG, const double *const Pos, const double *const PosOffset)
double	get_temp (double density, double ienergy, double helium, const struct UVBG *uvbg, double *ne_init)
double	GetNeutralFraction (double u_old, double rho, const struct UVBG *uvbg, double ne)
double	GetHeliumIonFraction (int ion, double u_old, double rho, const struct UVBG *uvbg, double ne_init)

Function Documentation

DoCooling()

double DoCooling	(	double	redshift,
		double	u_old,
		double	rho,
		double	dt,
		struct UVBG *	uvbg,
		double *	ne_guess,
		double	Z,
		double	MinEgySpec,
		int	isHeIIIionized
	)

Definition at line 71 of file cooling.c.

{
    if(!coolunits.CoolingOn) return 0;
 
    double u, du;
    double u_lower, u_upper;
    double LambdaNet;
    int iter = 0;
 
    rho *= coolunits.density_in_phys_cgs / PROTONMASS;  /* convert to (physical) protons/cm^3 */
    u_old *= coolunits.uu_in_cgs;
    MinEgySpec *= coolunits.uu_in_cgs;
    if(u_old < MinEgySpec)
        u_old = MinEgySpec;
    dt *= coolunits.tt_in_s;
 
    u = u_old;
    u_lower = u;
    u_upper = u;
 
    LambdaNet = get_lambdanet(rho, u, redshift, Z, uvbg, ne_guess, isHeIIIionized);
 
    /* bracketing */
    if(u - u_old - LambdaNet * dt < 0)  /* heating */
    {
        do
        {
            u_lower = u_upper;
            u_upper *= 1.1;
        } while(u_upper - u_old - get_lambdanet(rho, u_upper, redshift, Z, uvbg, ne_guess, isHeIIIionized) * dt < 0);
    }
    else
    {
        do {
            u_upper = u_lower;
            u_lower /= 1.1;
            /* This means that we don't need an initial bracket*/
            if(u_upper <= MinEgySpec) {
                break;
            }
        } while(u_lower - u_old - get_lambdanet(rho, u_lower, redshift, Z, uvbg, ne_guess, isHeIIIionized) * dt > 0);
    }
 
    do
    {
        u = 0.5 * (u_lower + u_upper);
        /* If we know that the new energy
         * is below the minimum gas internal energy, we are done here.*/
        if(u_upper <= MinEgySpec) {
            u = MinEgySpec;
            break;
        }
 
        LambdaNet = get_lambdanet(rho, u, redshift, Z, uvbg, ne_guess, isHeIIIionized);
 
        if(u - u_old - LambdaNet * dt > 0)
        {
            u_upper = u;
        }
        else
        {
            u_lower = u;
        }
 
        du = u_upper - u_lower;
 
        iter++;
 
        if(iter >= (MAXITER - 10))
            message(1, "u= %g\n", u);
    }
    while(fabs(du / u) > 1.0e-6 && iter < MAXITER);
 
    if(iter >= MAXITER)
    {
        endrun(10, "failed to converge in DoCooling()\n");
    }
 
    u /= coolunits.uu_in_cgs;   /*convert back to internal units */
 
    return u;
}

References cooling_units::CoolingOn, coolunits, cooling_units::density_in_phys_cgs, endrun(), get_lambdanet(), MAXITER, message(), PROTONMASS, cooling_units::tt_in_s, and cooling_units::uu_in_cgs.

Referenced by cooling_direct(), and test_DoCooling().

Here is the call graph for this function:

Here is the caller graph for this function:

get_global_UVBG()

struct UVBG get_global_UVBG

(

double

redshift

)

Definition at line 244 of file cooling_rates.c.

{
    struct UVBG GlobalUVBG = {0};
 
    if(!CoolingParams.PhotoIonizationOn)
        return GlobalUVBG;
 
    /* Set the homogeneous reionization redshift as the point when the UVB switches on.*/
    GlobalUVBG.zreion = pow(10,Gamma_log1z[NTreeCool - 1])-1;
 
    if(CoolingParams.UVRedshiftThreshold >= 0.)
        GlobalUVBG.zreion = CoolingParams.UVRedshiftThreshold;
 
    /* if a threshold is set, disable UV bg above that redshift */
    if(CoolingParams.UVRedshiftThreshold >= 0. && redshift > CoolingParams.UVRedshiftThreshold)
        return GlobalUVBG;
 
 
    GlobalUVBG.gJH0 = get_photo_rate(redshift, &Gamma_HI);
    GlobalUVBG.gJHe0 = get_photo_rate(redshift, &Gamma_HeI);
    GlobalUVBG.gJHep = get_photo_rate(redshift, &Gamma_HeII);
 
    GlobalUVBG.epsH0 = get_photo_rate(redshift, &Eps_HI);
    GlobalUVBG.epsHe0 = get_photo_rate(redshift, &Eps_HeI);
    /* During helium reionization we have a model for the inhomogeneous non-equilibrium heating.
     * To avoid double counting, remove the heating in the existing UVB*/
    if(during_helium_reionization(redshift))
        GlobalUVBG.epsHep = 0;
    else
        GlobalUVBG.epsHep = get_photo_rate(redshift, &Eps_HeII);
    GlobalUVBG.self_shield_dens = self_shield_dens(redshift, &GlobalUVBG);
    return GlobalUVBG;
}

References CoolingParams, cooling_params::fBar, UVBG::gJH0, GrayOpac, GrayOpac_ydata, GrayOpac_zz, and NGRAY.

Referenced by cooling_and_starformation(), get_helium_neutral_fraction_sfreff(), get_neutral_fraction_sfreff(), sfreff_on_eeqos(), test_DoCooling(), test_heatingcooling_rate(), test_rate_network(), and test_uvbg_loader().

Here is the caller graph for this function:

get_local_UVBG()

struct UVBG get_local_UVBG	(	double	redshift,
		const struct UVBG *const	GlobalUVBG,
		const double *const	Pos,
		const double *const	PosOffset
	)

Definition at line 91 of file cooling_uvfluc.c.

{
    if(!UVF.enabled) {
        /* directly use the TREECOOL table if UVF is disabled */
        return *GlobalUVBG;
    }
 
    struct UVBG uvbg = {0};
 
    uvbg.self_shield_dens = GlobalUVBG->self_shield_dens;
 
    double corrpos[3];
    int i;
    for(i = 0; i < 3; i++)
        corrpos[i] = Pos[i] - PosOffset[i];
    double zreion = interp_eval_periodic(&UVF.interp, corrpos, UVF.Table);
    if(zreion < redshift) {
        uvbg.zreion = zreion;
        return uvbg;
    }
    memcpy(&uvbg, GlobalUVBG, sizeof(struct UVBG));
    uvbg.zreion = zreion;
    return uvbg;
}

References CM_PER_MPC, endrun(), interp_init(), interp_init_dim(), message(), read_big_array(), UnitLength_in_cm, and UVF.

Referenced by cooling_direct(), get_helium_neutral_fraction_sfreff(), get_neutral_fraction_sfreff(), and get_sfr_eeqos().

Here is the call graph for this function:

Here is the caller graph for this function:

get_temp()

double get_temp	(	double	density,
		double	ienergy,
		double	helium,
		const struct UVBG *	uvbg,
		double *	ne_init
	)

Definition at line 1172 of file cooling_rates.c.

{
    double logt;
    double ne = get_equilib_ne(density, ienergy, helium, &logt, uvbg, *ne_init);
    double nh = density * (1 - helium);
    *ne_init = ne/nh;
    return get_temp_internal(ne/nh, ienergy, helium);
}

References density(), get_equilib_ne(), and get_temp_internal().

Referenced by test_rate_network().

Here is the call graph for this function:

Here is the caller graph for this function:

GetCoolingTime()

double GetCoolingTime	(	double	redshift,
		double	u_old,
		double	rho,
		struct UVBG *	uvbg,
		double *	ne_guess,
		double	Z
	)

Definition at line 157 of file cooling.c.

{
    if(!coolunits.CoolingOn) return 0;
 
    /* convert to physical cgs units */
    rho *= coolunits.density_in_phys_cgs / PROTONMASS;
    u_old *= coolunits.uu_in_cgs;
 
    /* Note: this does not include the long mean free path heating from helium reionization*/
    double LambdaNet = get_heatingcooling_rate(rho, u_old, 1 - HYDROGEN_MASSFRAC, redshift, Z, uvbg, ne_guess);
 
    if(LambdaNet >= 0)      /* ups, we have actually heating due to UV background */
        return 0;
 
    double coolingtime = u_old / (- LambdaNet);
 
    /*Convert back to internal units*/
    coolingtime /= coolunits.tt_in_s;
 
    return coolingtime;
}

References cooling_units::CoolingOn, coolunits, cooling_units::density_in_phys_cgs, get_heatingcooling_rate(), HYDROGEN_MASSFRAC, PROTONMASS, cooling_units::tt_in_s, and cooling_units::uu_in_cgs.

Referenced by cooling_relaxed(), get_egyeff(), init_cooling_and_star_formation(), and test_DoCooling().

Here is the call graph for this function:

Here is the caller graph for this function:

GetHeliumIonFraction()

double GetHeliumIonFraction	(	int	ion,
		double	u_old,
		double	rho,
		const struct UVBG *	uvbg,
		double	ne_init
	)

Definition at line 192 of file cooling.c.

{
    /* convert to physical cgs units */
    rho *= coolunits.density_in_phys_cgs / PROTONMASS;
    u_old *= coolunits.uu_in_cgs;
    double helium_ion = get_helium_ion_phys_cgs(ion, rho, u_old, 1 - HYDROGEN_MASSFRAC, uvbg, ne_init);
    return helium_ion;
}

References coolunits, cooling_units::density_in_phys_cgs, get_helium_ion_phys_cgs(), HYDROGEN_MASSFRAC, PROTONMASS, and cooling_units::uu_in_cgs.

Referenced by get_helium_neutral_fraction_sfreff().

Here is the call graph for this function:

Here is the caller graph for this function:

GetNeutralFraction()

double GetNeutralFraction	(	double	u_old,
		double	rho,
		const struct UVBG *	uvbg,
		double	ne
	)

Definition at line 181 of file cooling.c.

{
    /* convert to physical cgs units */
    rho *= coolunits.density_in_phys_cgs / PROTONMASS;
    u_old *= coolunits.uu_in_cgs;
    double nh0 = get_neutral_fraction_phys_cgs(rho, u_old, 1 - HYDROGEN_MASSFRAC, uvbg, &ne_init);
    return nh0;
}

References coolunits, cooling_units::density_in_phys_cgs, get_neutral_fraction_phys_cgs(), HYDROGEN_MASSFRAC, PROTONMASS, and cooling_units::uu_in_cgs.

Referenced by get_neutral_fraction_sfreff().

Here is the call graph for this function:

Here is the caller graph for this function:

init_cooling()

void init_cooling	(	const char *	TreeCoolFile,
		const char *	MetalCoolFile,
		char *	reion_hist_file,
		struct cooling_units	cu,
		Cosmology *	CP
	)

Definition at line 36 of file cooling.c.

{
    coolunits = cu;
    /*Initialize the cooling rates*/
    if(coolunits.CoolingOn)
        init_cooling_rates(TreeCoolFile, MetalCoolFile, CP);
    /* Initialize the helium reionization model*/
    init_qso_lightup(reion_hist_file);
}

References cooling_units::CoolingOn, coolunits, CP, init_cooling_rates(), and init_qso_lightup().

Referenced by init_cooling_and_star_formation(), and test_DoCooling().

Here is the call graph for this function:

Here is the caller graph for this function:

init_uvf_table()

void init_uvf_table	(	const char *	UVFluctuationFile,
		const int	UVFlucLen,
		const double	BoxSize,
		const double	UnitLength_in_cm
	)

Definition at line 91 of file cooling_uvfluc.c.

{
    if(strnlen(UVFluctuationFile, UVFlucLen) == 0) {
        UVF.enabled = 0;
        return;
    }
 
    /* Open and validate the UV fluctuation file*/
    BigFile bf;
    BigBlock bh;
    if(0 != big_file_mpi_open(&bf, UVFluctuationFile, MPI_COMM_WORLD)) {
        endrun(0, "Failed to open snapshot at %s:%s\n", UVFluctuationFile,
                    big_file_get_error_message());
    }
 
    if(0 != big_file_mpi_open_block(&bf, &bh, "Zreion_Table", MPI_COMM_WORLD)) {
        endrun(0, "Failed to create block at %s:%s\n", "Header",
                    big_file_get_error_message());
    }
    double TableBoxSize;
    double ReionRedshift;
    if ((0 != big_block_get_attr(&bh, "Nmesh", &UVF.Nside, "u8", 1)) ||
        (0 != big_block_get_attr(&bh, "BoxSize", &TableBoxSize, "f8", 1)) ||
        (0 != big_block_get_attr(&bh, "Redshift", &ReionRedshift, "f8", 1)) ||
        (0 != big_block_mpi_close(&bh, MPI_COMM_WORLD))) {
        endrun(0, "Failed to close block: %s\n",
                    big_file_get_error_message());
    }
    big_file_mpi_close(&bf, MPI_COMM_WORLD);
    double BoxMpc = BoxSize * UnitLength_in_cm / CM_PER_MPC;
    if(fabs(TableBoxSize - BoxMpc) > BoxMpc * 1e-5)
        endrun(0, "Wrong UV fluctuation file! %s is for box size %g Mpc/h, but current box is %g Mpc/h\n", UVFluctuationFile, TableBoxSize, BoxMpc);
 
    message(0, "Using NON-UNIFORM UV BG fluctuations from %s. Median reionization redshift is %g\n", UVFluctuationFile, ReionRedshift);
    UVF.enabled = 1;
 
    int size;
    UVF.Table = read_big_array(UVFluctuationFile, "Zreion_Table", &size);
 
    if(UVF.Nside * UVF.Nside * UVF.Nside != size)
        endrun(0, "Corrupt UV Fluctuation table: Nside = %ld, but table is %ld != %ld^3\n", UVF.Nside, size, UVF.Nside);
 
    int64_t dims[] = {UVF.Nside, UVF.Nside, UVF.Nside};
    interp_init(&UVF.interp, 3, dims);
    interp_init_dim(&UVF.interp, 0, 0, BoxSize);
    interp_init_dim(&UVF.interp, 1, 0, BoxSize);
    interp_init_dim(&UVF.interp, 2, 0, BoxSize);
 
    if(UVF.Table[0] < 0.01 || UVF.Table[0] > 100.0) {
        endrun(0, "UV Fluctuation out of range: %g\n", UVF.Table[0]);
    }
}

Referenced by init_cooling_and_star_formation().

Here is the caller graph for this function: