test_cosmology.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 <gsl/gsl_sf_hyperg.h>
#include <libgadget/physconst.h>
#include <libgadget/cosmology.h>
#include "stub.h"

Include dependency graph for test_cosmology.c:

Go to the source code of this file.

Functions
void	init_omega_nu (_omega_nu *const omnu, const double MNu[], const double a0, const double HubbleParam, const double tcmb0)
double	get_omega_nu (const _omega_nu *omnu, const double a)
double	get_omega_nu_nopart (const _omega_nu *omnu, const double a)
void	init_hybrid_nu (_hybrid_nu *const hybnu, const double mnu[], const double vcrit, const double light, const double nu_crit_time, const double kBtnu)
void	setup_cosmology (Cosmology *CP, double Omega0, double OmegaBaryon, double H0)
static double	radgrow (double aa, double omegar)
static double	growth (double aa, double omegam)
static void	test_cosmology (void **state)
int	main (void)

Function Documentation

get_omega_nu()

double get_omega_nu	(	const _omega_nu *const	omnu,
		const double	a
	)

Return the total matter density in neutrinos at scale factor a.

Definition at line 19 of file test_cosmology.c.

{
    return 0;
}

get_omega_nu_nopart()

double get_omega_nu_nopart	(	const _omega_nu *const	omnu,
		const double	a
	)

Return the total matter density in neutrinos at scale factor a , excluding active particles.

Definition at line 23 of file test_cosmology.c.

{
    return 0;
}

growth()

static double growth ( double  aa, double  omegam  )

inlinestatic

Definition at line 53 of file test_cosmology.c.

                                                      {
    double omegal = 1-omegam;
    return aa * gsl_sf_hyperg_2F1(1./3, 1, 11./6, -omegal/omegam*pow(aa,3));
}

Referenced by test_cosmology().

Here is the caller graph for this function:

init_hybrid_nu()

void init_hybrid_nu	(	_hybrid_nu *const	hybnu,
		const double	mnu[],
		const double	vcrit,
		const double	light,
		const double	nu_crit_time,
		const double	kBtnu
	)

Set up parameters for the hybrid neutrinos

Parameters

hybnu	initialised structure
mnu	array of neutrino masses in eV
vcrit	Critical velocity above which to treat neutrinos with particles. Note this is unperturbed velocity TODAY To get velocity at redshift z, multiply by (1+z)
light	speed of light in internal units
nu_crit_time	critical time to make neutrino particles live
kBtnu	Boltzmann constant times neutrino temperature. Dimensionful factor.

Definition at line 28 of file test_cosmology.c.

{ }

init_omega_nu()

void init_omega_nu	(	_omega_nu *const	omnu,
		const double	MNu[],
		const double	a0,
		const double	HubbleParam,
		const double	tcmb0
	)

Initialise the neutrino structure, do the time integration and allocate memory for the subclass rho_nu_single

Parameters

omnu	structure to initialise
MNu	array of neutrino masses in eV. Three entries.
a0	initial scale factor.
HubbleParam	Hubble parameter h0, eg, 0.7.
tcmb0	Redshift zero CMB temperature.

Definition at line 17 of file test_cosmology.c.

{}

main()

int main

(

void

)

Definition at line 96 of file test_cosmology.c.

               {
    const struct CMUnitTest tests[] = {
        cmocka_unit_test(test_cosmology),
    };
    return cmocka_run_group_tests_mpi(tests, NULL, NULL);
}

References test_cosmology().

Here is the call graph for this function:

radgrow()

static double radgrow ( double  aa, double  omegar  )

inlinestatic

Definition at line 48 of file test_cosmology.c.

                                                       {
    return omegar + 1.5 * 1. * aa;
}

Referenced by test_cosmology().

Here is the caller graph for this function:

setup_cosmology()

void setup_cosmology	(	Cosmology *	CP,
		double	Omega0,
		double	OmegaBaryon,
		double	H0
	)

Definition at line 30 of file test_cosmology.c.

{
    CP->CMBTemperature = 2.7255;
    CP->Omega0 = Omega0;
    CP->OmegaLambda = 1- CP->Omega0;
    CP->OmegaBaryon = OmegaBaryon;
    CP->HubbleParam = H0;
    CP->RadiationOn = 1;
    CP->Omega_fld = 0; /*Energy density of dark energy fluid at z=0*/
    CP->w0_fld = -1; /*Dark energy equation of state parameter*/
    CP->wa_fld = 0; /*Dark energy equation of state evolution parameter*/
    CP->Omega_ur = 0;
    CP->MNu[0] = CP->MNu[1] = CP->MNu[2] = 0;
    struct UnitSystem units = get_unitsystem(3.085678e21, 1.989e43, 1e5);
    /*Do the main cosmology initialisation*/
    init_cosmology(CP,0.01, units);
}

References Cosmology::CMBTemperature, CP, get_unitsystem(), Cosmology::HubbleParam, init_cosmology(), Cosmology::MNu, Cosmology::Omega0, Cosmology::Omega_fld, Cosmology::Omega_ur, Cosmology::OmegaBaryon, Cosmology::OmegaLambda, Cosmology::RadiationOn, Cosmology::w0_fld, and Cosmology::wa_fld.

Referenced by test_cosmology().

Here is the call graph for this function:

Here is the caller graph for this function:

test_cosmology()

static void test_cosmology ( void **  state )

static

Definition at line 58 of file test_cosmology.c.

{
    Cosmology CP = {0};
    //Check that we get the right scalings for total matter domination.
    //Cosmology(double HubbleParam, double Omega, double OmegaLambda, double MNu, int Hierarchy, bool NoRadiation)
    setup_cosmology(&CP, 1., 0.0455, 0.7);
    /*Check some of the setup worked*/
    assert_true(CP.OmegaK < 1e-5);
    assert_true(fabs(CP.OmegaG/5.045e-5 - 1) < 2e-3);
    /*Check the hubble function is sane*/
    CP.RadiationOn = 0;
    assert_true(fabs(hubble_function(&CP, 1) - CP.Hubble) < 1e-5);
    CP.RadiationOn = 1;
    assert_true(fabs(hubble_function(&CP, 1) - CP.Hubble* sqrt(1+CP.OmegaG)) < 1e-7);
 
    assert_true(fabs(CP.Hubble - 0.1) < 1e-6);
    assert_true((hubble_function(&CP, 1) - CP.Hubble) < 1e-5);
    assert_true(fabs(hubble_function(&CP, 0.1) - hubble_function(&CP, 1)/pow(0.1,3/2.)) < 1e-2);
    assert_true(fabs(GrowthFactor(&CP, 0.5,1.)/0.5 -1) < 2e-4);
    //Check that the velocity correction d ln D1/d lna is constant
    assert_true(fabs(1.0 - F_Omega(&CP, 1.5)) < 1e-1);
    assert_true(fabs(1.0 - F_Omega(&CP, 2)) < 1e-2);
    //Check radiation against exact solution from gr-qc/0504089
    assert_true(fabs(1/GrowthFactor(&CP, 0.05,1.) - radgrow(1., CP.OmegaG)/radgrow(0.05, CP.OmegaG))< 1e-3);
    assert_true(fabs(GrowthFactor(&CP, 0.01,0.001) - radgrow(0.01, CP.OmegaG)/radgrow(0.001, CP.OmegaG))< 1e-3);
 
    //Check against exact solutions from gr-qc/0504089: No radiation!
    //Note that the GSL hyperg needs the last argument to be < 1
    double omegam = 0.5;
    setup_cosmology(&CP, omegam, 0.0455, 0.7);
    CP.RadiationOn = 0;
    //Check growth factor during matter domination
    assert_true(fabs(1/GrowthFactor(&CP, 0.5, 1.) - growth(1., omegam)/growth(0.5, omegam)) < 1e-3);
    assert_true(fabs(GrowthFactor(&CP, 0.3, 0.15) - growth(0.3, omegam)/growth(0.15, omegam)) < 1e-3);
    assert_true(fabs(1/GrowthFactor(&CP, 0.01, 1.) - growth(1, omegam)/growth(0.01, omegam)) < 1e-3);
    assert_true(fabs(0.01*log(GrowthFactor(&CP, 0.01+1e-5,0.01-1e-5))/2e-5 -  F_Omega(&CP, 0.01)) < 1e-3);
}

References CP, F_Omega(), growth(), GrowthFactor(), Cosmology::Hubble, hubble_function(), Cosmology::OmegaG, Cosmology::OmegaK, radgrow(), Cosmology::RadiationOn, and setup_cosmology().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function: