test_forcetree.c File Reference

#include <stdarg.h>
#include <stddef.h>
#include <setjmp.h>
#include <cmocka.h>
#include <math.h>
#include <mpi.h>
#include <string.h>
#include <stdio.h>
#include <time.h>
#include <omp.h>
#include <gsl/gsl_rng.h>
#include <libgadget/forcetree.h>
#include <libgadget/partmanager.h>
#include <libgadget/domain.h>
#include "stub.h"

Include dependency graph for test_forcetree.c:

Go to the source code of this file.

Classes
struct	forcetree_testdata

Macros
#define	NODECACHE_SIZE   100

Functions
void	dump_snapshot ()
double	walltime_measure_full (char *name, char *file, int line)
static int	order_by_type_and_key (const void *a, const void *b)
static int	check_moments (const ForceTree *tb, const int numpart, const int nrealnode)
static int	check_tree (const ForceTree *tb, const int nnodes, const int numpart)
static void	do_tree_test (const int numpart, ForceTree tb, DomainDecomp *ddecomp)
static void	test_rebuild_flat (void **state)
static void	test_rebuild_close (void **state)
void	do_random_test (gsl_rng *r, const int numpart, const ForceTree tb, DomainDecomp *ddecomp)
static void	test_rebuild_random (void **state)
void	trivial_domain (DomainDecomp *ddecomp)
static int	setup_tree (void **state)
static int	teardown_tree (void **state)
int	main (void)

Variables
struct part_manager_type	PartManager [1]
struct slots_manager_type	SlotsManager [1]

Macro Definition Documentation

NODECACHE_SIZE

#define NODECACHE_SIZE   100

Definition at line 63 of file test_forcetree.c.

Function Documentation

check_moments()

static int check_moments ( const ForceTree *  tb, const int  numpart, const int  nrealnode  )

static

Definition at line 67 of file test_forcetree.c.

{
    double * oldmass = malloc(sizeof(double) * tb->numnodes);
    int i;
 
    for(i=tb->firstnode; i < tb->numnodes + tb->firstnode; i ++) {
        oldmass[i - tb->firstnode] = tb->Nodes[i].mom.mass;
    }
 
    for(i=0; i<numpart; i++)
    {
        int fnode = force_get_father(i, tb);
        /*Subtract mass so that nothing is left.*/
        assert_true(fnode >= tb->firstnode && fnode < tb->lastnode);
        while(fnode > 0) {
            tb->Nodes[fnode].mom.mass -= P[i].Mass;
            fnode = tb->Nodes[fnode].father;
            /*Validate father*/
            assert_true((fnode >= tb->firstnode && fnode < tb->lastnode) || fnode == -1);
        }
    }
    int node = tb->firstnode;
    int counter = 0;
    int sibcntr = 0;
    while(node >= 0) {
        /* Assert a real node*/
        assert_true(node >= -1 && node < tb->lastnode && node >= tb->firstnode);
        struct NODE * nop = &tb->Nodes[node];
 
        /*Check sibling*/
        assert_true(tb->Nodes[node].sibling >= -1 && tb->Nodes[node].sibling < tb->lastnode);
        int sib = tb->Nodes[node].sibling;
        int sfather = force_get_father(sib, tb);
        int father = force_get_father(node, tb);
        /* Our sibling should either be a true sibling, with the same father,
            * or should be the child of one of our ancestors*/
        if(sfather != father && sib != -1) {
            int ances = father;
            while(ances >= 0) {
                assert_true(ances >= tb->firstnode);
                ances = force_get_father(ances, tb);
                if(ances == sfather)
                    break;
            }
            assert_int_equal(ances, sfather);
/*                 printf("node %d ances %d sib %d next %d father %d sfather %d\n",node, ances, sib, nop->s.suns[0], father, sfather); */
        }
        else if(sib == -1)
            sibcntr++;
 
        if(!(tb->Nodes[node].mom.mass < 0.5 && tb->Nodes[node].mom.mass > -0.5)) {
            printf("node %d (%d) mass %g / %g TL %d DLM %d ITL %d\n",
                node, node - tb->firstnode, tb->Nodes[node].mom.mass, oldmass[node - tb->firstnode],
                tb->Nodes[node].f.TopLevel,
                tb->Nodes[node].f.DependsOnLocalMass,
                tb->Nodes[node].f.InternalTopLevel
            );
            /* something is wrong show the particles */
            if(tb->Nodes[node].f.ChildType == PARTICLE_NODE_TYPE)
                for(i = 0; i < nop->s.noccupied; i++) {
                    int nn = nop->s.suns[i];
                    printf("particles P[%d], Mass=%g\n", nn, P[nn].Mass);
                }
        }
        assert_true(tb->Nodes[node].mom.mass < 0.5 && tb->Nodes[node].mom.mass > -0.5);
        /*Check center of mass moments*/
        for(i=0; i<3; i++)
            assert_true(tb->Nodes[node].mom.cofm[i] <= PartManager->BoxSize && tb->Nodes[node].mom.cofm[i] >= 0);
        counter++;
 
        if(nop->f.ChildType == PARTICLE_NODE_TYPE)
            node = nop->sibling;
        else
            node = nop->s.suns[0];
    }
//     message(5, "count %d real %d\n", counter, nrealnode);
    assert_true(counter <= nrealnode);
    assert_true(sibcntr < counter/100);
 
    free(oldmass);
    return nrealnode;
}

References part_manager_type::BoxSize, NODE::ChildType, NODE::cofm, NODE::DependsOnLocalMass, NODE::f, NODE::father, ForceTree::firstnode, force_get_father(), NODE::InternalTopLevel, ForceTree::lastnode, NODE::mass, NODE::mom, NodeChild::noccupied, ForceTree::Nodes, ForceTree::numnodes, P, PARTICLE_NODE_TYPE, PartManager, NODE::s, NODE::sibling, NodeChild::suns, and NODE::TopLevel.

Referenced by do_tree_test().

Here is the call graph for this function:

Here is the caller graph for this function:

check_tree()

static int check_tree ( const ForceTree *  tb, const int  nnodes, const int  numpart  )

static

Definition at line 153 of file test_forcetree.c.

{
    const int firstnode = tb->firstnode;
    int tot_empty = 0, nrealnode = 0, sevens = 0;
    int i;
    for(i=firstnode; i<nnodes+firstnode; i++)
    {
        struct NODE * pNode = &(tb->Nodes[i]);
        /*Just reserved free space with nothing in it*/
        if(pNode->father < -1.5)
            continue;
 
        int j;
        /* Full of particles*/
        if(pNode->s.noccupied < 1<<16) {
            tot_empty += NMAXCHILD - pNode->s.noccupied;
            if(pNode->s.noccupied == 0)
                sevens++;
            for(j=0; j<pNode->s.noccupied; j++) {
                int child = pNode->s.suns[j];
                assert_true(child >= 0);
                assert_true(child < firstnode);
                P[child].PI += 1;
                assert_int_equal(force_get_father(child, tb), i);
            }
        }
        /* Node is full of other nodes*/
        else {
            for(j=0; j<8; j++) {
                /*Check children*/
                int child = pNode->s.suns[j];
                assert_true(child < firstnode+nnodes);
                assert_true(child >= firstnode);
                assert_true(fabs(tb->Nodes[child].len/pNode->len - 0.5) < 1e-4);
                int k;
                for(k=0; k<3; k++) {
                    if(j & (1<<k))
                        assert_true(tb->Nodes[child].center[k] > pNode->center[k]);
                    else
                        assert_true(tb->Nodes[child].center[k] <= pNode->center[k]);
                }
            }
        }
        nrealnode++;
    }
 
    for(i=0; i<numpart; i++)
    {
        assert_int_equal(P[i].PI, 1);
    }
    printf("Tree filling factor: %g on %d nodes (wasted: %d empty: %d)\n", tot_empty/(8.*nrealnode), nrealnode, nnodes - nrealnode, sevens);
    return nrealnode - sevens;
}

References NODE::center, NODE::father, ForceTree::firstnode, force_get_father(), NODE::len, NMAXCHILD, NodeChild::noccupied, ForceTree::Nodes, P, NODE::s, and NodeChild::suns.

Referenced by do_tree_test().

Here is the call graph for this function:

Here is the caller graph for this function:

do_random_test()

void do_random_test	(	gsl_rng *	r,
		const int	numpart,
		const ForceTree	tb,
		DomainDecomp *	ddecomp
	)

Definition at line 300 of file test_forcetree.c.

{
    /* Create a regular grid of particles, 8x8x8, all of type 1,
     * in a box 8 kpc across.*/
    int i;
    for(i=0; i<numpart/4; i++) {
        P[i].Type = 1;
        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 = 1;
        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 = 1;
        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_tree_test(numpart, tb, ddecomp);
}

References part_manager_type::BoxSize, forcetree_testdata::ddecomp, do_tree_test(), P, PartManager, and forcetree_testdata::r.

Referenced by test_rebuild_random().

Here is the call graph for this function:

Here is the caller graph for this function:

do_tree_test()

static void do_tree_test ( const int  numpart, ForceTree  tb, DomainDecomp *  ddecomp  )

static

Definition at line 207 of file test_forcetree.c.

{
    /*Sort by peano key so this is more realistic*/
    int i;
    #pragma omp parallel for
    for(i=0; i<numpart; i++) {
        P[i].Key = PEANO(P[i].Pos, PartManager->BoxSize);
        P[i].Mass = 1;
        P[i].PI = 0;
        P[i].IsGarbage = 0;
    }
    qsort(P, numpart, sizeof(struct particle_data), order_by_type_and_key);
    int maxnode = tb.lastnode - tb.firstnode;
    PartManager->MaxPart = numpart;
    PartManager->NumPart = numpart;
    assert_true(tb.Nodes != NULL);
    /*So we know which nodes we have initialised*/
    for(i=0; i< maxnode; i++)
        tb.Nodes_base[i].father = -2;
    /*Time creating the nodes*/
    double start, end;
    start = MPI_Wtime();
    int nodes = force_tree_create_nodes(tb, numpart, ALLMASK, ddecomp, 0);
    tb.numnodes = nodes;
    assert_true(nodes < maxnode);
    end = MPI_Wtime();
    double ms = (end - start)*1000;
    printf("Number of nodes used: %d. Built tree in %.3g ms\n", nodes,ms);
    int nrealnode = check_tree(&tb, nodes, numpart);
    /* now compute the multipole moments recursively */
    start = MPI_Wtime();
    force_update_node_parallel(&tb, ddecomp);
    end = MPI_Wtime();
    ms = (end - start)*1000;
    printf("Updated moments in %.3g ms. Total mass: %g\n", ms, tb.Nodes[tb.firstnode].mom.mass);
    assert_true(fabs(tb.Nodes[tb.firstnode].mom.mass - numpart) < 0.5);
    check_moments(&tb, numpart, nrealnode);
}

References ALLMASK, part_manager_type::BoxSize, check_moments(), check_tree(), NODE::father, ForceTree::firstnode, force_tree_create_nodes(), force_update_node_parallel(), ForceTree::lastnode, NODE::mass, part_manager_type::MaxPart, NODE::mom, ForceTree::Nodes, ForceTree::Nodes_base, ForceTree::numnodes, part_manager_type::NumPart, order_by_type_and_key(), P, PartManager, and PEANO().

Referenced by do_random_test(), test_rebuild_close(), and test_rebuild_flat().

Here is the call graph for this function:

Here is the caller graph for this function:

dump_snapshot()

void dump_snapshot

(

)

Definition at line 37 of file test_forcetree.c.

{ }

main()

int main

(

void

)

Definition at line 397 of file test_forcetree.c.

               {
    const struct CMUnitTest tests[] = {
        cmocka_unit_test(test_rebuild_flat),
        cmocka_unit_test(test_rebuild_close),
        cmocka_unit_test(test_rebuild_random),
    };
    return cmocka_run_group_tests_mpi(tests, setup_tree, teardown_tree);
}

References setup_tree(), teardown_tree(), test_rebuild_close(), test_rebuild_flat(), and test_rebuild_random().

Here is the call graph for this function:

order_by_type_and_key()

static int order_by_type_and_key ( const void *  a, const void *  b  )

static

Definition at line 46 of file test_forcetree.c.

{
    const struct particle_data * pa  = (const struct particle_data *) a;
    const struct particle_data * pb  = (const struct particle_data *) b;
 
    if(pa->Type < pb->Type)
        return -1;
    if(pa->Type > pb->Type)
        return +1;
    if(pa->Key < pb->Key)
        return -1;
    if(pa->Key > pb->Key)
        return +1;
 
    return 0;
}

References particle_data::Key, and particle_data::Type.

Referenced by do_tree_test().

Here is the caller graph for this function:

setup_tree()

static int setup_tree ( void **  state )

static

Definition at line 371 of file test_forcetree.c.

                                    {
    /*Set up the important parts of the All structure.*/
    /*Particles should not be outside this*/
    memset(PartManager, 0, sizeof(PartManager[0]));
    memset(SlotsManager, 0, sizeof(SlotsManager[0]));
    PartManager->BoxSize = 8;
    init_forcetree_params(2);
    /*Set up the top-level domain grid*/
    struct forcetree_testdata *data = malloc(sizeof(struct forcetree_testdata));
    trivial_domain(&data->ddecomp);
    data->r = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(data->r, 0);
    *state = (void *) data;
    return 0;
}

References part_manager_type::BoxSize, forcetree_testdata::ddecomp, init_forcetree_params(), PartManager, forcetree_testdata::r, SlotsManager, and trivial_domain().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

teardown_tree()

static int teardown_tree ( void **  state )

static

Definition at line 387 of file test_forcetree.c.

                                       {
    struct forcetree_testdata * data = (struct forcetree_testdata * ) *state;
    free(data->ddecomp.TopNodes);
    free(data->ddecomp.TopLeaves);
    free(data->ddecomp.Tasks);
    free(data->r);
    free(data);
    return 0;
}

References forcetree_testdata::ddecomp, forcetree_testdata::r, DomainDecomp::Tasks, DomainDecomp::TopLeaves, and DomainDecomp::TopNodes.

Referenced by main().

Here is the caller graph for this function:

test_rebuild_close()

static void test_rebuild_close ( void **  state )

static

Definition at line 276 of file test_forcetree.c.

                                              {
    /*Set up the particle data*/
    int ncbrt = 128;
    int numpart = ncbrt*ncbrt*ncbrt;
    double close = 5000;
    P = malloc(numpart*sizeof(struct particle_data));
    /* Create particles clustered in one place, all of type 1.*/
    int i;
    #pragma omp parallel for
    for(i=0; i<numpart; i++) {
        P[i].Type = 1;
        P[i].Pos[0] = 4. + (i/ncbrt/ncbrt)/close;
        P[i].Pos[1] = 4. + ((i/ncbrt) % ncbrt) /close;
        P[i].Pos[2] = 4. + (i % ncbrt)/close;
    }
    struct forcetree_testdata * data = * (struct forcetree_testdata **) state;
    DomainDecomp ddecomp = data->ddecomp;
    ddecomp.TopLeaves[0].topnode = numpart;
    ForceTree tb = force_treeallocate(numpart, numpart, &ddecomp);
    do_tree_test(numpart, tb, &ddecomp);
    force_tree_free(&tb);
    free(P);
}

References forcetree_testdata::ddecomp, do_tree_test(), force_tree_free(), force_treeallocate(), P, DomainDecomp::TopLeaves, and topleaf_data::topnode.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

test_rebuild_flat()

static void test_rebuild_flat ( void **  state )

static

Definition at line 246 of file test_forcetree.c.

                                             {
    /*Set up the particle data*/
    int ncbrt = 128;
    int numpart = ncbrt*ncbrt*ncbrt;
    P = malloc(numpart*sizeof(struct particle_data));
    /* 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 = 1;
        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);
    }
    /*Allocate tree*/
    /*Base pointer*/
    struct forcetree_testdata * data = * (struct forcetree_testdata **) state;
    DomainDecomp ddecomp = data->ddecomp;
    ddecomp.TopLeaves[0].topnode = numpart;
    ForceTree tb = force_treeallocate(numpart, numpart, &ddecomp);
    /* So unused memory has Father < 0*/
    for(i = tb.firstnode; i < tb.lastnode; i++)
        tb.Nodes[i].father = -10;
 
    do_tree_test(numpart, tb, &ddecomp);
    force_tree_free(&tb);
    free(P);
}

References part_manager_type::BoxSize, forcetree_testdata::ddecomp, do_tree_test(), NODE::father, ForceTree::firstnode, force_tree_free(), force_treeallocate(), ForceTree::lastnode, ForceTree::Nodes, P, PartManager, DomainDecomp::TopLeaves, and topleaf_data::topnode.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

test_rebuild_random()

static void test_rebuild_random ( void **  state )

static

Definition at line 326 of file test_forcetree.c.

                                               {
    /*Set up the particle data*/
    int ncbrt = 64;
    struct forcetree_testdata * data = * (struct forcetree_testdata **) state;
    DomainDecomp ddecomp = data->ddecomp;
    gsl_rng * r = (gsl_rng *) data->r;
    int numpart = ncbrt*ncbrt*ncbrt;
    /*Allocate tree*/
    /*Base pointer*/
    ddecomp.TopLeaves[0].topnode = numpart;
    ForceTree tb = force_treeallocate(numpart, numpart, &ddecomp);
    assert_true(tb.Nodes != NULL);
    P = malloc(numpart*sizeof(struct particle_data));
    int i;
    for(i=0; i<2; i++) {
        do_random_test(r, numpart, tb, &ddecomp);
    }
    force_tree_free(&tb);
    free(P);
}

References forcetree_testdata::ddecomp, do_random_test(), force_tree_free(), force_treeallocate(), P, forcetree_testdata::r, DomainDecomp::TopLeaves, and topleaf_data::topnode.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

trivial_domain()

void trivial_domain

(

DomainDecomp *

ddecomp

)

Definition at line 348 of file test_forcetree.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 = malloc(sizeof(struct topnode_data));
    ddecomp->TopNodes[0].Daughter = -1;
    ddecomp->TopNodes[0].Leaf = 0;
    ddecomp->TopLeaves = malloc(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 = malloc(sizeof(struct task_data));
    ddecomp->Tasks[0].StartLeaf = 0;
    ddecomp->Tasks[0].EndLeaf = 1;
}

References BITS_PER_DIMENSION, topnode_data::Daughter, forcetree_testdata::ddecomp, DomainDecomp::domain_allocated_flag, task_data::EndLeaf, topnode_data::Leaf, 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_tree().

Here is the caller graph for this function:

walltime_measure_full()

double walltime_measure_full	(	char *	name,
		char *	file,
		int	line
	)

Definition at line 39 of file test_forcetree.c.

                                                                 {
    return MPI_Wtime();
}

Variable Documentation

PartManager

struct part_manager_type PartManager[1]

Definition at line 1 of file test_forcetree.c.

Referenced by check_moments(), do_random_test(), do_tree_test(), setup_tree(), and test_rebuild_flat().

SlotsManager

struct slots_manager_type SlotsManager[1]

Definition at line 1 of file test_forcetree.c.

Referenced by setup_tree().