test_cooling_rates.c File Reference

#include <stdarg.h>
#include <stddef.h>
#include <setjmp.h>
#include <cmocka.h>
#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <libgadget/physconst.h>
#include <libgadget/cooling_rates.h>
#include <libgadget/config.h>
#include "stub.h"

Include dependency graph for test_cooling_rates.c:

Go to the source code of this file.

Macros
#define	NEXACT   20

Functions
int	during_helium_reionization (double redshift)
double	recomb_alphaHp (double temp)
static struct cooling_params	get_test_coolpar (void)
static void	test_recomb_rates (void **state)
static void	test_uvbg_loader (void **state)
static void	test_rate_network (void **state)
static void	test_heatingcooling_rate (void **state)
int	main (void)

Variables
static const double	f92g2 [NEXACT] = {5.758e-11, 2.909e-11, 1.440e-11, 6.971e-12,3.282e-12, 1.489e-12, 6.43e-13, 2.588e-13, 9.456e-14, 3.069e-14, 8.793e-15, 2.245e-15, 5.190e-16, 1.107e-16, 2.221e-17, 4.267e-18, 7.960e-19, 1.457e-19,2.636e-20, 4.737e-21}
static const double	f92n1 [NEXACT] = {9.258e-12, 5.206e-12, 2.927e-12, 1.646e-12, 9.246e-13, 5.184e-13, 2.890e-13, 1.582e-13, 8.255e-14, 3.882e-14, 1.545e-14, 5.058e-15, 1.383e-15, 3.276e-16, 7.006e-17, 1.398e-17, 2.665e-18, 4.940e-19, 9.001e-20, 1.623e-20}
static const double	tt [NEXACT] = {3.16227766e+00, 1.0e+01, 3.16227766e+01, 1.0e+02, 3.16227766e+02, 1.00e+03, 3.16227766e+03, 1.e+04, 3.16227766e+04, 1.e+05, 3.16227766e+05, 1.e+06, 3.16227766e+06, 1.0e+07, 3.16227766e+07, 1.0e+08, 3.16227766e+08, 1.0e+09, 3.16227766e+09, 1.0e+10}

Macro Definition Documentation

NEXACT

#define NEXACT   20

Definition at line 25 of file test_cooling_rates.c.

Function Documentation

during_helium_reionization()

int during_helium_reionization

(

double

redshift

)

Definition at line 18 of file test_cooling_rates.c.

{
    return 0;
}

get_test_coolpar()

static struct cooling_params get_test_coolpar ( void  )

static

Definition at line 32 of file test_cooling_rates.c.

{
    static struct cooling_params coolpar;
    coolpar.CMBTemperature = 2.7255;
    coolpar.PhotoIonizeFactor = 1;
    coolpar.SelfShieldingOn = 1;
    coolpar.fBar = 0.17;
    coolpar.PhotoIonizationOn = 1;
    coolpar.recomb = Verner96;
    coolpar.cooling = Sherwood;
    coolpar.UVRedshiftThreshold = -1;
    coolpar.MinGasTemp = 100;
    coolpar.HeliumHeatOn = 0;
    coolpar.HydrogenHeatAmp = 0;
    return coolpar;
}

Referenced by test_heatingcooling_rate(), test_rate_network(), test_recomb_rates(), and test_uvbg_loader().

Here is the caller graph for this function:

main()

int main

(

void

)

Definition at line 321 of file test_cooling_rates.c.

               {
    const struct CMUnitTest tests[] = {
        cmocka_unit_test(test_recomb_rates),
        cmocka_unit_test(test_rate_network),
        cmocka_unit_test(test_heatingcooling_rate),
        cmocka_unit_test(test_uvbg_loader)
    };
    return cmocka_run_group_tests_mpi(tests, NULL, NULL);
}

References test_heatingcooling_rate(), test_rate_network(), test_recomb_rates(), and test_uvbg_loader().

Here is the call graph for this function:

recomb_alphaHp()

double recomb_alphaHp

(

double

temp

)

Definition at line 347 of file cooling_rates.c.

{
    switch(CoolingParams.recomb)
    {
        case Cen92:
            return 8.4e-11 / sqrt(temp) / pow(temp/1000, 0.2) / (1+ pow(temp/1e6, 0.7));
        case Verner96:
            /*The V&F 96 fitting formula is accurate to < 1% in the worst case.
            See line 1 of V&F96 table 1.*/
            return _Verner96Fit(temp, 7.982e-11, 0.748, 3.148, 7.036e+05);
        case Badnell06:
            /*Slightly different fit coefficients*/
            return _Verner96Fit(temp, 8.318e-11, 0.7472, 2.965, 7.001e5);
        default:
            endrun(3, "Recombination rate not supported\n");
    }
}

References _Verner96Fit(), Badnell06, Cen92, CoolingParams, endrun(), cooling_params::recomb, and Verner96.

Referenced by cool_RecombHp(), init_cooling_rates(), nH0_internal(), recomb_alphaHepp(), and test_recomb_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

test_heatingcooling_rate()

static void test_heatingcooling_rate ( void **  state )

static

Definition at line 174 of file test_cooling_rates.c.

{
    struct cooling_params coolpar = get_test_coolpar();
    coolpar.recomb = Cen92;
    coolpar.cooling = KWH92;
    coolpar.SelfShieldingOn = 0;
 
    const char * TreeCool = GADGET_TESTDATA_ROOT "/examples/TREECOOL_ep_2018p";
    const char * MetalCool = "";
 
    /*unit system*/
    double HubbleParam = 0.697;
    double UnitDensity_in_cgs = 6.76991e-22;
    double UnitTime_in_s = 3.08568e+16;
    double UnitMass_in_g = 1.989e+43;
    double UnitLength_in_cm = 3.08568e+21;
    double UnitEnergy_in_cgs = UnitMass_in_g  * pow(UnitLength_in_cm, 2) / pow(UnitTime_in_s, 2);
    Cosmology CP = {0};
    CP.OmegaCDM = 0.3;
    CP.OmegaBaryon = coolpar.fBar * CP.OmegaCDM;
    CP.HubbleParam = HubbleParam;
 
    set_coolpar(coolpar);
    init_cooling_rates(TreeCool, MetalCool, &CP);
 
    struct cooling_units coolunits;
    coolunits.CoolingOn = 1;
    coolunits.density_in_phys_cgs = UnitDensity_in_cgs * HubbleParam * HubbleParam;
    coolunits.uu_in_cgs = UnitEnergy_in_cgs / UnitMass_in_g;
    coolunits.tt_in_s = UnitTime_in_s / HubbleParam;
 
    /*Default values from sfr_eff.c. Some dependence on HubbleParam, so don't change it.*/
    double egyhot = 2104.92;
    egyhot *= coolunits.uu_in_cgs;
 
    /* convert to physical cgs units */
    double dens = 0.027755;
    dens *= coolunits.density_in_phys_cgs/PROTONMASS;
    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. */
    double LambdaNet = get_heatingcooling_rate(dens, egyhot, 1 - HYDROGEN_MASSFRAC, 0, 0, &uvbg, &ne);
 
    double tcool = egyhot / (- LambdaNet);
 
    /*Convert back to internal units*/
    tcool /= coolunits.tt_in_s;
 
    //message(1, "tcool = %g LambdaNet = %g ne=%g\n", tcool, LambdaNet, ne);
    /* This differs by 0.13% from the old cooling code number,
     * apparently just because of rounding errors. The excitation cooling
     * numbers from Cen are not accurate to better than 1% anyway, so don't worry about it*/
    assert_true(fabs(tcool / 4.68906e-06 - 1) < 1e-3);
 
    /*Now check that we get the desired cooling rate with a UVB*/
    uvbg = get_global_UVBG(0);
 
    assert_true(uvbg.epsHep > 0);
    assert_true(uvbg.gJHe0 > 0);
 
    dens /= 100;
    LambdaNet = get_heatingcooling_rate(dens, egyhot/10., 1 - HYDROGEN_MASSFRAC, 0, 0, &uvbg, &ne);
    //message(1, "LambdaNet = %g, uvbg=%g\n", LambdaNet, uvbg.epsHep);
    assert_true(fabs(LambdaNet / (-0.0410059) - 1) < 1e-3);
 
    LambdaNet = get_heatingcooling_rate(dens/2.5, egyhot/10., 1 - HYDROGEN_MASSFRAC, 0, 0, &uvbg, &ne);
    assert_true(LambdaNet > 0);
    /*Check self-shielding affects the cooling rates*/
    coolpar.SelfShieldingOn = 1;
    set_coolpar(coolpar);
    init_cooling_rates(TreeCool, MetalCool, &CP);
    LambdaNet = get_heatingcooling_rate(dens*1.5, egyhot/10., 1 - HYDROGEN_MASSFRAC, 0, 0, &uvbg, &ne);
    //message(1, "LambdaNet = %g, uvbg=%g\n", LambdaNet, uvbg.epsHep);
    assert_false(LambdaNet > 0);
    assert_true(fabs(LambdaNet/ (-1.64834) - 1) < 1e-3);
}

References Cen92, cooling_params::cooling, cooling_units::CoolingOn, coolunits, CP, cooling_units::density_in_phys_cgs, UVBG::epsHep, cooling_params::fBar, get_global_UVBG(), get_heatingcooling_rate(), get_test_coolpar(), UVBG::gJHe0, Cosmology::HubbleParam, HYDROGEN_MASSFRAC, init_cooling_rates(), KWH92, MetalCool, Cosmology::OmegaBaryon, Cosmology::OmegaCDM, PROTONMASS, cooling_params::recomb, cooling_params::SelfShieldingOn, set_coolpar(), cooling_units::tt_in_s, UnitLength_in_cm, and cooling_units::uu_in_cgs.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

test_rate_network()

static void test_rate_network ( void **  state )

static

Definition at line 119 of file test_cooling_rates.c.

{
    struct cooling_params coolpar = get_test_coolpar();
    const char * TreeCool = GADGET_TESTDATA_ROOT "/examples/TREECOOL_ep_2018p";
    const char * MetalCool = "";
 
    set_coolpar(coolpar);
    Cosmology CP = {0};
    CP.OmegaCDM = 0.3;
    CP.OmegaBaryon = coolpar.fBar * CP.OmegaCDM;
    CP.HubbleParam = 0.7;
 
    init_cooling_rates(TreeCool, MetalCool, &CP);
 
    struct UVBG uvbg = get_global_UVBG(2);
 
 
    //Complete ionisation at low density
    double logt;
    assert_true( fabs(get_equilib_ne(1e-6, 200.*1e10, 0.24, &logt, &uvbg, 1) / (1e-6*0.76) - (1 + 2* 0.24/(1-0.24)/4)) < 3e-5);
    assert_true( fabs(get_equilib_ne(1e-6, 200.*1e10, 0.12, &logt, &uvbg, 1) / (1e-6*0.88) - (1 + 2* 0.12/(1-0.12)/4)) < 3e-5);
    assert_true( fabs(get_equilib_ne(1e-5, 200.*1e10, 0.24, &logt, &uvbg, 1) / (1e-5*0.76) - (1 + 2* 0.24/(1-0.24)/4)) < 3e-4);
    assert_true( fabs(get_equilib_ne(1e-4, 200.*1e10, 0.24, &logt, &uvbg, 1) / (1e-4*0.76) - (1 + 2* 0.24/(1-0.24)/4)) < 2e-3);
 
    double ne = 1.;
    double temp = get_temp(1e-4, 200.*1e10,0.24, &uvbg, &ne);
    assert_true(9500 < temp);
    assert_true(temp < 9510);
    //Roughly prop to internal energy when ionised
    assert_true(fabs(get_temp(1e-4, 400.*1e10,0.24, &uvbg, &ne) / get_temp(1e-4, 200.*1e10,0.24, &uvbg, &ne) - 2.) < 1e-3);
 
    assert_true(fabs(get_temp(1, 200.*1e10,0.24, &uvbg, &ne) - 14700) < 200);
 
    //Neutral fraction prop to density.
    double dens[3] = {1e-4, 1e-5, 1e-6};
    int i;
    for(i = 0; i < 3; i++) {
        assert_true(fabs(get_neutral_fraction_phys_cgs(dens[i], 200.*1e10,0.24, &uvbg, &ne) / dens[i] - 0.3113) < 1e-3);
    }
    //Neutral (self-shielded) at high density:
    assert_true(get_neutral_fraction_phys_cgs(1, 100.,0.24, &uvbg, &ne) > 0.95);
    assert_true(0.75 > get_neutral_fraction_phys_cgs(0.1, 100.*1e10,0.24, &uvbg, &ne));
    assert_true(get_neutral_fraction_phys_cgs(0.1, 100.*1e10,0.24, &uvbg, &ne) > 0.735);
 
    //Check self-shielding is working.
    coolpar.SelfShieldingOn = 0;
    set_coolpar(coolpar);
    init_cooling_rates(TreeCool, MetalCool, &CP);
 
    assert_true( get_neutral_fraction_phys_cgs(1, 100.*1e10,0.24, &uvbg, &ne) < 0.25);
    assert_true( get_neutral_fraction_phys_cgs(0.1, 100.*1e10,0.24, &uvbg, &ne) <0.05);
}

References CP, cooling_params::fBar, get_equilib_ne(), get_global_UVBG(), get_neutral_fraction_phys_cgs(), get_temp(), get_test_coolpar(), Cosmology::HubbleParam, init_cooling_rates(), MetalCool, Cosmology::OmegaBaryon, Cosmology::OmegaCDM, cooling_params::SelfShieldingOn, and set_coolpar().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

test_recomb_rates()

static void test_recomb_rates ( void **  state )

static

Definition at line 52 of file test_cooling_rates.c.

{
    struct cooling_params coolpar = get_test_coolpar();
    int i;
    const char * TreeCool = GADGET_TESTDATA_ROOT "/examples/TREECOOL_ep_2018p";
    const char * MetalCool = "";
 
    Cosmology CP = {0};
    CP.OmegaCDM = 0.3;
    CP.OmegaBaryon = coolpar.fBar * CP.OmegaCDM;
    CP.HubbleParam = 0.7;
 
    set_coolpar(coolpar);
    init_cooling_rates(TreeCool, MetalCool, &CP);
    for(i=0; i< NEXACT; i++) {
        assert_true(fabs(recomb_alphaHp(tt[i])/(f92g2[i]+f92n1[i])-1.) < 1e-2);
    }
 
    coolpar.recomb = Cen92;
    set_coolpar(coolpar);
    init_cooling_rates(TreeCool, MetalCool, &CP);
    /*Cen rates are not very accurate.*/
    for(i=4; i< 12; i++) {
        assert_true(fabs(recomb_alphaHp(tt[i])/(f92g2[i]+f92n1[i])-1.) < 0.5);
    }
}

References Cen92, CP, f92g2, f92n1, cooling_params::fBar, get_test_coolpar(), Cosmology::HubbleParam, init_cooling_rates(), MetalCool, NEXACT, Cosmology::OmegaBaryon, Cosmology::OmegaCDM, cooling_params::recomb, recomb_alphaHp(), set_coolpar(), and tt.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

test_uvbg_loader()

static void test_uvbg_loader ( void **  state )

static

Definition at line 80 of file test_cooling_rates.c.

{
    struct cooling_params coolpar = get_test_coolpar();
    coolpar.SelfShieldingOn = 1;
    const char * TreeCool = GADGET_TESTDATA_ROOT "/examples/TREECOOL_ep_2018p";
    const char * MetalCool = "";
    set_coolpar(coolpar);
    Cosmology CP = {0};
    CP.OmegaCDM = 0.3;
    CP.OmegaBaryon = coolpar.fBar * CP.OmegaCDM;
    CP.HubbleParam = 0.7;
    init_cooling_rates(TreeCool, MetalCool, &CP);
    /*Test sensible rates at high redshift*/
    struct UVBG uvbg = get_global_UVBG(16);
    assert_true(uvbg.epsH0 == 0);
    assert_true(uvbg.self_shield_dens > 1e8);
    assert_true(uvbg.gJH0 == 0);
    /*Test at zero redshift*/
    uvbg = get_global_UVBG(0);
    assert_true(fabs(uvbg.epsH0/3.65296e-25 -1) < 1e-5);
    assert_true(fabs(uvbg.epsHe0/3.98942e-25 -1) < 1e-5);
    assert_true(fabs(uvbg.epsHep/3.33253e-26 -1) < 1e-5);
    assert_true(fabs(uvbg.gJH0/6.06e-14 -1) < 1e-5);
    assert_true(fabs(uvbg.gJHe0/3.03e-14 -1) < 1e-5);
    assert_true(fabs(uvbg.gJHep/1.1e-15 -1) < 1e-5);
    assert_true(fabs(uvbg.self_shield_dens/0.0010114161149989826 - 1) < 1e-5);
    /*Test at intermediate redshift*/
    uvbg = get_global_UVBG(3.);
    //message(0, "uvbg %g %g %g %g %g %g %g\n", uvbg.gJH0, uvbg.gJHe0, uvbg.gJHep,  uvbg.epsH0, uvbg.epsHe0, uvbg.epsHep, uvbg.self_shield_dens);
    assert_true(fabs(uvbg.epsH0/5.96570906168362e-24 -1) < 1e-5);
    assert_true(fabs(uvbg.epsHe0/4.466976578202419e-24 -1) < 1e-5);
    assert_true(fabs(uvbg.epsHep/2.758535690259892e-26 -1) < 1e-5);
    assert_true(fabs(uvbg.gJH0/1.0549960730284017e-12 -1) < 1e-5);
    assert_true(fabs(uvbg.gJHe0/4.759025257653999e-13 -1) < 1e-5);
    assert_true(fabs(uvbg.gJHep/2.270599708640625e-16 -1) < 1e-5);
    assert_true(fabs(uvbg.self_shield_dens/0.007691709693529007 - 1) < 1e-5);
}

References CP, UVBG::epsH0, UVBG::epsHe0, UVBG::epsHep, cooling_params::fBar, get_global_UVBG(), get_test_coolpar(), UVBG::gJH0, UVBG::gJHe0, UVBG::gJHep, Cosmology::HubbleParam, init_cooling_rates(), MetalCool, Cosmology::OmegaBaryon, Cosmology::OmegaCDM, UVBG::self_shield_dens, cooling_params::SelfShieldingOn, and set_coolpar().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

f92g2

const double f92g2[NEXACT] = {5.758e-11, 2.909e-11, 1.440e-11, 6.971e-12,3.282e-12, 1.489e-12, 6.43e-13, 2.588e-13, 9.456e-14, 3.069e-14, 8.793e-15, 2.245e-15, 5.190e-16, 1.107e-16, 2.221e-17, 4.267e-18, 7.960e-19, 1.457e-19,2.636e-20, 4.737e-21}

static

Definition at line 29 of file test_cooling_rates.c.

Referenced by test_recomb_rates().

f92n1

const double f92n1[NEXACT] = {9.258e-12, 5.206e-12, 2.927e-12, 1.646e-12, 9.246e-13, 5.184e-13, 2.890e-13, 1.582e-13, 8.255e-14, 3.882e-14, 1.545e-14, 5.058e-15, 1.383e-15, 3.276e-16, 7.006e-17, 1.398e-17, 2.665e-18, 4.940e-19, 9.001e-20, 1.623e-20}

static

Definition at line 31 of file test_cooling_rates.c.

Referenced by test_recomb_rates().

tt

const double tt[NEXACT] = {3.16227766e+00, 1.0e+01, 3.16227766e+01, 1.0e+02, 3.16227766e+02, 1.00e+03, 3.16227766e+03, 1.e+04, 3.16227766e+04, 1.e+05, 3.16227766e+05, 1.e+06, 3.16227766e+06, 1.0e+07, 3.16227766e+07, 1.0e+08, 3.16227766e+08, 1.0e+09, 3.16227766e+09, 1.0e+10}

static

Definition at line 32 of file test_cooling_rates.c.

Referenced by init_cooling_rates(), and test_recomb_rates().