test_density.c File Reference

#include <stdarg.h>
#include <stddef.h>
#include <setjmp.h>
#include <cmocka.h>
#include <math.h>
#include <mpi.h>
#include <stdio.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <libgadget/partmanager.h>
#include <libgadget/walltime.h>
#include <libgadget/slotsmanager.h>
#include <libgadget/utils/mymalloc.h>
#include <libgadget/density.h>
#include <libgadget/domain.h>
#include <libgadget/forcetree.h>
#include <libgadget/timestep.h>
#include <libgadget/gravity.h>
#include "stub.h"

Include dependency graph for test_density.c:

Go to the source code of this file.

Classes
struct	density_testdata

Functions
static void	set_init_hsml (ForceTree *Tree)
static void	check_densities (double MinGasHsml)
static void	do_density_test (void **state, const int numpart, double expectedhsml, double hsmlerr)
static void	test_density_flat (void **state)
static void	test_density_close (void **state)
void	do_random_test (void **state, gsl_rng *r, const int numpart)
static void	test_density_random (void **state)
void	trivial_domain (DomainDecomp *ddecomp)
static int	teardown_density (void **state)
static int	setup_density (void **state)
int	main (void)

Variables
static struct ClockTable	CT

Function Documentation

check_densities()

static void check_densities ( double  MinGasHsml )

static

Definition at line 60 of file test_density.c.

{
    int i;
    double maxHsml=P[0].Hsml, minHsml= P[0].Hsml;
    #pragma omp parallel for reduction(min:minHsml) reduction(max:maxHsml)
    for(i=0; i<PartManager->NumPart; i++) {
        assert_true(isfinite(P[i].Hsml));
        assert_true(isfinite(SPHP(i).Density));
        assert_true(SPHP(i).Density > 0);
        if(P[i].Hsml < minHsml)
            minHsml = P[i].Hsml;
        if(P[i].Hsml > maxHsml)
            maxHsml = P[i].Hsml;
    }
    assert_true(isfinite(minHsml));
    assert_true(minHsml >= MinGasHsml);
    assert_true(maxHsml <= PartManager->BoxSize);
 
}

References part_manager_type::NumPart, P, PartManager, and SPHP.

Referenced by do_density_test().

Here is the caller graph for this function:

do_density_test()

static void do_density_test ( void **  state, const int  numpart, double  expectedhsml, double  hsmlerr  )

static

Definition at line 80 of file test_density.c.

{
    int i, npbh=0;
    #pragma omp parallel for reduction(+: npbh)
    for(i=0; i<numpart; i++) {
        int j;
        P[i].Key = PEANO(P[i].Pos, PartManager->BoxSize);
        P[i].Mass = 1;
        P[i].TimeBin = 0;
        P[i].Ti_drift = 0;
        for(j=0; j<3; j++)
            P[i].Vel[j] = 1.5;
        if(P[i].Type == 0) {
            SPHP(i).Entropy = 1;
            SPHP(i).DtEntropy = 0;
            SPHP(i).Density = 1;
        }
        if(P[i].Type == 5)
            npbh++;
    }
 
    SlotsManager->info[0].size = numpart-npbh;
    SlotsManager->info[5].size = npbh;
    PartManager->NumPart = numpart;
    ActiveParticles act = {0};
    act.NumActiveParticle = numpart;
    act.ActiveParticle = NULL;
    struct density_testdata * data = * (struct density_testdata **) state;
    DomainDecomp ddecomp = data->ddecomp;
    ddecomp.TopLeaves[0].topnode = PartManager->MaxPart;
 
    ForceTree tree = {0};
    force_tree_rebuild(&tree, &ddecomp, 0, 1, NULL);
    set_init_hsml(&tree);
    /* Rebuild without moments to check it works*/
    force_tree_rebuild(&tree, &ddecomp, 0, 0, NULL);
    /*Time doing the density finding*/
    double start, end;
    start = MPI_Wtime();
    /*Find the density*/
    DriftKickTimes kick = {0};
    Cosmology CP = {0};
    CP.CMBTemperature = 2.7255;
    CP.Omega0 = 0.3;
    CP.OmegaLambda = 1- CP.Omega0;
    CP.OmegaBaryon = 0.045;
    CP.HubbleParam = 0.7;
    CP.RadiationOn = 0;
    CP.w0_fld = -1; /*Dark energy equation of state parameter*/
    /*Should be 0.1*/
    struct UnitSystem units = get_unitsystem(3.085678e21, 1.989e43, 1e5);
    init_cosmology(&CP,0.01, units);
 
    density(&act, 1, 0, 0, 0, kick, &CP, &data->sph_pred, NULL, &tree);
    end = MPI_Wtime();
    double ms = (end - start)*1000;
    message(0, "Found densities in %.3g ms\n", ms);
    check_densities(data->dp.MinGasHsmlFractional);
 
    double avghsml = 0;
    #pragma omp parallel for reduction(+:avghsml)
    for(i=0; i<numpart; i++) {
        avghsml += P[i].Hsml;
    }
    message(0, "Average Hsml: %g Expected %g +- %g\n",avghsml/numpart, expectedhsml, hsmlerr);
    assert_true(fabs(avghsml/numpart - expectedhsml) < hsmlerr);
    /* Make MaxNumNgbDeviation smaller and check we get a consistent result.*/
    double * Hsml = mymalloc2("Hsml", numpart * sizeof(double));
    #pragma omp parallel for
    for(i=0; i<numpart; i++) {
        Hsml[i] = P[i].Hsml;
    }
    data->dp.MaxNumNgbDeviation = 0.5;
    set_densitypar(data->dp);
 
    start = MPI_Wtime();
    /*Find the density*/
    density(&act, 1, 0, 0, 0, kick, &CP, &data->sph_pred, NULL, &tree);
    end = MPI_Wtime();
    ms = (end - start)*1000;
    message(0, "Found 1 dev densities in %.3g ms\n", ms);
    double diff = 0;
    double DesNumNgb = GetNumNgb(GetDensityKernelType());
    /* Free tree before checks so that we still recover if checks fail*/
    force_tree_free(&tree);
 
    #pragma omp parallel for reduction(max:diff)
    for(i=0; i<numpart; i++) {
        assert_true(fabs(Hsml[i]/P[i].Hsml-1) < data->dp.MaxNumNgbDeviation / DesNumNgb);
        if(fabs(Hsml[i] - P[i].Hsml) > diff)
            diff = fabs(Hsml[i] - P[i].Hsml);
    }
    message(0, "Max diff between Hsml: %g\n",diff);
    myfree(Hsml);
 
    check_densities(data->dp.MinGasHsmlFractional);
}

References ActiveParticles::ActiveParticle, part_manager_type::BoxSize, check_densities(), Cosmology::CMBTemperature, CP, density_testdata::ddecomp, density(), density_testdata::dp, force_tree_free(), force_tree_rebuild(), get_unitsystem(), GetDensityKernelType(), GetNumNgb(), Cosmology::HubbleParam, slots_manager_type::info, init_cosmology(), density_params::MaxNumNgbDeviation, part_manager_type::MaxPart, message(), density_params::MinGasHsmlFractional, myfree, mymalloc2, ActiveParticles::NumActiveParticle, part_manager_type::NumPart, Cosmology::Omega0, Cosmology::OmegaBaryon, Cosmology::OmegaLambda, P, PartManager, PEANO(), Cosmology::RadiationOn, set_densitypar(), set_init_hsml(), slot_info::size, SlotsManager, density_testdata::sph_pred, SPHP, DomainDecomp::TopLeaves, topleaf_data::topnode, and Cosmology::w0_fld.

Referenced by do_random_test(), test_density_close(), and test_density_flat().

Here is the call graph for this function:

Here is the caller graph for this function:

do_random_test()

void do_random_test	(	void **	state,
		gsl_rng *	r,
		const int	numpart
	)

Definition at line 230 of file test_density.c.

{
    /* Create a randomly space set of particles, 8x8x8, all of type 0. */
    int i;
    for(i=0; i<numpart/4; i++) {
        P[i].Type = 0;
        P[i].PI = i;
        P[i].Hsml = PartManager->BoxSize/cbrt(numpart);
 
        int j;
        for(j=0; j<3; j++)
            P[i].Pos[j] = PartManager->BoxSize * gsl_rng_uniform(r);
    }
    for(i=numpart/4; i<3*numpart/4; i++) {
        P[i].Type = 0;
        P[i].PI = i;
        P[i].Hsml = PartManager->BoxSize/cbrt(numpart);
        int j;
        for(j=0; j<3; j++)
            P[i].Pos[j] = PartManager->BoxSize/2 + PartManager->BoxSize/8 * exp(pow(gsl_rng_uniform(r)-0.5,2));
    }
    for(i=3*numpart/4; i<numpart; i++) {
        P[i].Type = 0;
        P[i].PI = i;
        P[i].Hsml = PartManager->BoxSize/cbrt(numpart);
        int j;
        for(j=0; j<3; j++)
            P[i].Pos[j] = PartManager->BoxSize*0.1 + PartManager->BoxSize/32 * exp(pow(gsl_rng_uniform(r)-0.5,2));
    }
    do_density_test(state, numpart, 0.187515, 1e-3);
}

References part_manager_type::BoxSize, do_density_test(), P, and PartManager.

Referenced by test_density_random().

Here is the call graph for this function:

Here is the caller graph for this function:

main()

int main

(

void

)

Definition at line 357 of file test_density.c.

               {
    const struct CMUnitTest tests[] = {
        cmocka_unit_test(test_density_flat),
        cmocka_unit_test(test_density_close),
        cmocka_unit_test(test_density_random),
    };
    return cmocka_run_group_tests_mpi(tests, setup_density, teardown_density);
}

set_init_hsml()

static void set_init_hsml ( ForceTree *  Tree )

static

Definition at line 34 of file test_density.c.

{
    int i;
    #pragma omp parallel for
    for(i = 0; i < PartManager->NumPart; i++)
    {
        int no = force_get_father(i, Tree);
 
        double DesNumNgb = GetNumNgb(GetDensityKernelType());
        while(10 * DesNumNgb * P[i].Mass > Tree->Nodes[no].mom.mass)
        {
            int p = force_get_father(no, Tree);
 
            if(p < 0)
                break;
 
            no = p;
        }
 
        P[i].Hsml =
            pow(3.0 / (4 * M_PI) * DesNumNgb * P[i].Mass / (Tree->Nodes[no].mom.mass),
                    1.0 / 3) * Tree->Nodes[no].len;
    }
}

References force_get_father(), GetDensityKernelType(), GetNumNgb(), NODE::len, NODE::mass, NODE::mom, ForceTree::Nodes, part_manager_type::NumPart, P, and PartManager.

Referenced by do_density_test().

Here is the call graph for this function:

Here is the caller graph for this function:

setup_density()

static int setup_density ( void **  state )

static

Definition at line 314 of file test_density.c.

                                       {
    /* Needed so the integer timeline works*/
    setup_sync_points(0.01, 0.1, 0.0, 0);
 
    /*Reserve space for the slots*/
    slots_init(0.01 * PartManager->MaxPart, SlotsManager);
    slots_set_enabled(0, sizeof(struct sph_particle_data), SlotsManager);
    slots_set_enabled(5, sizeof(struct sph_particle_data), SlotsManager);
    int64_t atleast[6] = {0};
    atleast[0] = pow(32,3);
    atleast[5] = 2;
    int64_t maxpart = 0;
    int i;
    for(i = 0; i < 6; i++)
        maxpart+=atleast[i];
    const double BoxSize = 8;
    particle_alloc_memory(PartManager, BoxSize, maxpart);
    slots_reserve(1, atleast, SlotsManager);
    walltime_init(&CT);
    init_forcetree_params(2);
    struct density_testdata *data = mymalloc("data", sizeof(struct density_testdata));
    data->sph_pred = slots_allocate_sph_pred_data(maxpart);
    /*Set up the top-level domain grid*/
    trivial_domain(&data->ddecomp);
    data->dp.DensityResolutionEta = 1.;
    data->dp.BlackHoleNgbFactor = 2;
    data->dp.MaxNumNgbDeviation = 2;
    data->dp.DensityKernelType = DENSITY_KERNEL_CUBIC_SPLINE;
    data->dp.MinGasHsmlFractional = 0.006;
    struct gravshort_tree_params tree_params = {0};
    tree_params.FractionalGravitySoftening = 1;
    set_gravshort_treepar(tree_params);
 
    gravshort_set_softenings(1);
    data->dp.BlackHoleMaxAccretionRadius = 99999.;
 
    set_densitypar(data->dp);
    data->r = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(data->r, 0);
    *state = (void *) data;
    return 0;
}

References density_params::BlackHoleMaxAccretionRadius, density_params::BlackHoleNgbFactor, CT, density_testdata::ddecomp, DENSITY_KERNEL_CUBIC_SPLINE, density_params::DensityKernelType, density_params::DensityResolutionEta, density_testdata::dp, gravshort_tree_params::FractionalGravitySoftening, gravshort_set_softenings(), init_forcetree_params(), density_params::MaxNumNgbDeviation, part_manager_type::MaxPart, density_params::MinGasHsmlFractional, mymalloc, particle_alloc_memory(), PartManager, density_testdata::r, set_densitypar(), set_gravshort_treepar(), setup_sync_points(), slots_allocate_sph_pred_data(), slots_init(), slots_reserve(), slots_set_enabled(), SlotsManager, density_testdata::sph_pred, trivial_domain(), and walltime_init().

Here is the call graph for this function:

teardown_density()

static int teardown_density ( void **  state )

static

Definition at line 301 of file test_density.c.

                                          {
    struct density_testdata * data = (struct density_testdata * ) *state;
    slots_free_sph_pred_data(&data->sph_pred);
    myfree(data->ddecomp.Tasks);
    myfree(data->ddecomp.TopLeaves);
    myfree(data->ddecomp.TopNodes);
    free(data->r);
    myfree(data);
    return 0;
}

test_density_close()

static void test_density_close ( void **  state )

static

Definition at line 197 of file test_density.c.

                                              {
    /*Set up the particle data*/
    int ncbrt = 32;
    int numpart = ncbrt*ncbrt*ncbrt;
    double close = 500.;
    int i;
    /* A few particles scattered about the place so the tree is not sparse*/
    #pragma omp parallel for
    for(i=0; i<numpart/4; i++) {
        P[i].Type = 0;
        P[i].PI = i;
        P[i].Hsml = 4*PartManager->BoxSize/cbrt(numpart/8);
        P[i].Pos[0] = (PartManager->BoxSize/ncbrt) * (i/(ncbrt/2.)/(ncbrt/2.));
        P[i].Pos[1] = (PartManager->BoxSize/ncbrt) * ((i*2/ncbrt) % (ncbrt/2));
        P[i].Pos[2] = (PartManager->BoxSize/ncbrt) * (i % (ncbrt/2));
    }
 
    /* Create particles clustered in one place, all of type 0.*/
    #pragma omp parallel for
    for(i=numpart/4; i<numpart; i++) {
        P[i].Type = 0;
        P[i].PI = i;
        P[i].Hsml = 2*ncbrt/close;
        P[i].Pos[0] = 4.1 + (i/ncbrt/ncbrt)/close;
        P[i].Pos[1] = 4.1 + ((i/ncbrt) % ncbrt) /close;
        P[i].Pos[2] = 4.1 + (i % ncbrt)/close;
    }
    P[numpart-1].Type = 5;
    P[numpart-1].PI = 0;
 
    do_density_test(state, numpart, 0.131726, 1e-4);
}

References part_manager_type::BoxSize, do_density_test(), P, and PartManager.

Here is the call graph for this function:

test_density_flat()

static void test_density_flat ( void **  state )

static

Definition at line 178 of file test_density.c.

                                             {
    /*Set up the particle data*/
    int ncbrt = 32;
    int numpart = ncbrt*ncbrt*ncbrt;
    /* Create a regular grid of particles, 8x8x8, all of type 1,
     * in a box 8 kpc across.*/
    int i;
    #pragma omp parallel for
    for(i=0; i<numpart; i++) {
        P[i].Type = 0;
        P[i].PI = i;
        P[i].Hsml = 1.5*PartManager->BoxSize/cbrt(numpart);
        P[i].Pos[0] = (PartManager->BoxSize/ncbrt) * (i/ncbrt/ncbrt);
        P[i].Pos[1] = (PartManager->BoxSize/ncbrt) * ((i/ncbrt) % ncbrt);
        P[i].Pos[2] = (PartManager->BoxSize/ncbrt) * (i % ncbrt);
    }
    do_density_test(state, numpart, 0.501747, 1e-4);
}

References part_manager_type::BoxSize, do_density_test(), P, and PartManager.

Here is the call graph for this function:

test_density_random()

static void test_density_random ( void **  state )

static

Definition at line 262 of file test_density.c.

                                               {
    /*Set up the particle data*/
    int ncbrt = 32;
    struct density_testdata * data = * (struct density_testdata **) state;
    gsl_rng * r = (gsl_rng *) data->r;
    int numpart = ncbrt*ncbrt*ncbrt;
    /*Allocate tree*/
    /*Base pointer*/
    int i;
    for(i=0; i<2; i++) {
        do_random_test(state, r, numpart);
    }
}

References do_random_test(), and density_testdata::r.

Here is the call graph for this function:

trivial_domain()

void trivial_domain

(

DomainDecomp *

ddecomp

)

Definition at line 278 of file test_density.c.

{
    /* The whole tree goes into one topnode.
     * Set up just enough of the TopNode structure that
     * domain_get_topleaf works*/
    ddecomp->domain_allocated_flag = 1;
    ddecomp->NTopNodes = 1;
    ddecomp->NTopLeaves = 1;
    ddecomp->TopNodes = mymalloc("topnode", sizeof(struct topnode_data));
    ddecomp->TopNodes[0].Daughter = -1;
    ddecomp->TopNodes[0].Leaf = 0;
    ddecomp->TopLeaves = mymalloc("topleaf",sizeof(struct topleaf_data));
    ddecomp->TopLeaves[0].Task = 0;
    ddecomp->TopLeaves[0].topnode = 0;
    /*These are not used*/
    ddecomp->TopNodes[0].StartKey = 0;
    ddecomp->TopNodes[0].Shift = BITS_PER_DIMENSION * 3;
    /*To tell the code we are in serial*/
    ddecomp->Tasks = mymalloc("task",sizeof(struct task_data));
    ddecomp->Tasks[0].StartLeaf = 0;
    ddecomp->Tasks[0].EndLeaf = 1;
}

References BITS_PER_DIMENSION, topnode_data::Daughter, density_testdata::ddecomp, DomainDecomp::domain_allocated_flag, task_data::EndLeaf, topnode_data::Leaf, mymalloc, DomainDecomp::NTopLeaves, DomainDecomp::NTopNodes, topnode_data::Shift, topnode_data::StartKey, task_data::StartLeaf, topleaf_data::Task, DomainDecomp::Tasks, DomainDecomp::TopLeaves, topleaf_data::topnode, and DomainDecomp::TopNodes.

Referenced by setup_density().

Here is the caller graph for this function:

Variable Documentation

CT

struct ClockTable CT

static

Definition at line 301 of file test_density.c.

Referenced by setup_density().