test_gravity.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 <string.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <omp.h>
#include "stub.h"
#include <libgadget/utils/mymalloc.h>
#include <libgadget/utils/system.h>
#include <libgadget/utils/endrun.h>
#include <libgadget/partmanager.h>
#include <libgadget/walltime.h>
#include <libgadget/domain.h>
#include <libgadget/forcetree.h>
#include <libgadget/gravity.h>
#include <libgadget/petapm.h>
#include <libgadget/timestep.h>
#include <libgadget/physconst.h>

Include dependency graph for test_gravity.c:

Go to the source code of this file.

Classes
struct	forcetree_testdata

Functions
static void	grav_force (const int this, const int other, const double *offset, double *accns)
void	check_accns (double *meanerr_tot, double *maxerr_tot, double *PairAccn, double meanacc)
static void	find_means (double *meangrav, double *suppmean, double *suppaccns)
static void	force_direct (double *accn)
static int	check_against_force_direct (double ErrTolForceAcc)
static void	do_force_test (int Nmesh, double Asmth, double ErrTolForceAcc, int direct)
static void	test_force_flat (void **state)
static void	test_force_close (void **state)
void	do_random_test (gsl_rng *r, const int numpart)
static void	test_force_random (void **state)
static int	setup_tree (void **state)
static int	teardown_tree (void **state)
int	main (void)

Variables
static struct ClockTable	CT
static const double	G = 43.0071

Function Documentation

check_accns()

void check_accns	(	double *	meanerr_tot,
		double *	maxerr_tot,
		double *	PairAccn,
		double	meanacc
	)

Definition at line 73 of file test_gravity.c.

{
    double meanerr=0, maxerr=-1;
    int64_t i;
    /* This checks that the short-range force accuracy is being correctly estimated.*/
    #pragma omp parallel for reduction(+: meanerr) reduction(max: maxerr)
    for(i = 0; i < PartManager->NumPart; i++)
    {
        int k;
        for(k=0; k<3; k++) {
            double err = fabs((PairAccn[3*i+k] - (P[i].GravPM[k] + P[i].GravAccel[k]))/meanacc);
            meanerr += err;
            if(maxerr < err)
                maxerr = err;
        }
    }
    MPI_Allreduce(&meanerr, meanerr_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&maxerr, maxerr_tot, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
    int64_t tot_npart;
    MPI_Allreduce(&PartManager->NumPart, &tot_npart, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
 
    *meanerr_tot/= (tot_npart*3.);
}

References GravPM, MPI_INT64, part_manager_type::NumPart, P, and PartManager.

Referenced by check_against_force_direct().

Here is the caller graph for this function:

check_against_force_direct()

static int check_against_force_direct ( double  ErrTolForceAcc )

static

Definition at line 146 of file test_gravity.c.

{
    double * accn = (double *) mymalloc("accelerations", 3*sizeof(double) * PartManager->NumPart);
    force_direct(accn);
    double meanerr=0, maxerr=-1, meanacc=0, meanforce=0;
    find_means(&meanacc, &meanforce, accn);
    check_accns(&meanerr, &maxerr, accn, meanacc);
    myfree(accn);
    message(0, "Mean rel err is: %g max rel err is %g, meanacc %g mean grav force %g\n", meanerr, maxerr, meanacc, meanforce);
    /*Make some statements about the force error*/
    assert_true(maxerr < 3*ErrTolForceAcc);
    assert_true(meanerr < 0.8*ErrTolForceAcc);
 
    return 0;
}

References check_accns(), find_means(), force_direct(), message(), myfree, mymalloc, part_manager_type::NumPart, and PartManager.

Referenced by do_force_test().

Here is the call graph for this function:

Here is the caller graph for this function:

do_force_test()

static void do_force_test ( int  Nmesh, double  Asmth, double  ErrTolForceAcc, int  direct  )

static

Definition at line 162 of file test_gravity.c.

{
    /*Sort by peano key so this is more realistic*/
    int i;
    #pragma omp parallel for
    for(i=0; i<PartManager->NumPart; i++) {
        P[i].Type = 1;
        P[i].Key = PEANO(P[i].Pos, PartManager->BoxSize);
        P[i].Mass = 1;
        P[i].ID = i;
        P[i].TimeBin = 0;
        P[i].IsGarbage = 0;
    }
 
    ActiveParticles act = {0};
    act.NumActiveParticle = PartManager->NumPart;
 
    DomainDecomp ddecomp = {0};
    domain_decompose_full(&ddecomp);
 
    PetaPM pm = {0};
    gravpm_init_periodic(&pm, PartManager->BoxSize, Asmth, Nmesh, G);
    ForceTree Tree = {0};
    force_tree_rebuild(&Tree, &ddecomp, 1, 1, NULL);
    gravshort_fill_ntab(SHORTRANGE_FORCE_WINDOW_TYPE_EXACT, Asmth);
    /* Setup cosmology*/
    Cosmology CP ={0};
    CP.MNu[0] = CP.MNu[1] = CP.MNu[2] = 0;
    CP.OmegaCDM = 0.3;
    CP.CMBTemperature = 2.72;
    CP.HubbleParam = 0.7;
    CP.Omega0 = 0.3;
    CP.OmegaBaryon = 0.045;
    CP.OmegaLambda = 0.7;
    struct UnitSystem units = get_unitsystem(3.085678e21, 1.989e43, 1e5);
    init_cosmology(&CP, 0.01, units);
 
    gravpm_force(&pm, &Tree, &CP, 0.1, CM_PER_MPC/1000., ".", 0.01, 2);
    force_tree_rebuild(&Tree, &ddecomp, 1, 1, NULL);
    const double rho0 = CP.Omega0 * 3 * CP.Hubble * CP.Hubble / (8 * M_PI * G);
 
    /* Barnes-Hut on first iteration*/
    struct gravshort_tree_params treeacc = {0};
    treeacc.BHOpeningAngle = 0.175;
    treeacc.TreeUseBH = 1;
    treeacc.Rcut = 7;
    treeacc.ErrTolForceAcc = ErrTolForceAcc;
    treeacc.AdaptiveSoftening = 0;
    treeacc.FractionalGravitySoftening = 1./30.;
 
    set_gravshort_treepar(treeacc);
    gravshort_set_softenings(PartManager->BoxSize / cbrt(PartManager->NumPart));
 
    /* Twice so the opening angle is consistent*/
    grav_short_tree(&act, &pm, &Tree, rho0, 0, 2);
    grav_short_tree(&act, &pm, &Tree, rho0, 0, 2);
 
    force_tree_free(&Tree);
    petapm_destroy(&pm);
    domain_free(&ddecomp);
    if(direct)
        check_against_force_direct(ErrTolForceAcc);
}

References gravshort_tree_params::AdaptiveSoftening, gravshort_tree_params::BHOpeningAngle, part_manager_type::BoxSize, check_against_force_direct(), CM_PER_MPC, Cosmology::CMBTemperature, CP, domain_decompose_full(), domain_free(), gravshort_tree_params::ErrTolForceAcc, force_tree_free(), force_tree_rebuild(), gravshort_tree_params::FractionalGravitySoftening, G, get_unitsystem(), grav_short_tree(), gravpm_force(), gravpm_init_periodic(), gravshort_fill_ntab(), gravshort_set_softenings(), Cosmology::Hubble, Cosmology::HubbleParam, init_cosmology(), Cosmology::MNu, ActiveParticles::NumActiveParticle, part_manager_type::NumPart, Cosmology::Omega0, Cosmology::OmegaBaryon, Cosmology::OmegaCDM, Cosmology::OmegaLambda, P, PartManager, PEANO(), petapm_destroy(), gravshort_tree_params::Rcut, set_gravshort_treepar(), SHORTRANGE_FORCE_WINDOW_TYPE_EXACT, and gravshort_tree_params::TreeUseBH.

Referenced by do_random_test(), test_force_close(), and test_force_flat().

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
	)

Definition at line 291 of file test_gravity.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++) {
        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++) {
        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++) {
        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));
    }
    PartManager->NumPart = numpart;
    PartManager->MaxPart = numpart;
    do_force_test(48, 1.5, 0.002, 1);
}

References part_manager_type::BoxSize, do_force_test(), part_manager_type::MaxPart, part_manager_type::NumPart, P, and PartManager.

Referenced by test_force_random().

Here is the call graph for this function:

Here is the caller graph for this function:

find_means()

static void find_means ( double *  meangrav, double *  suppmean, double *  suppaccns  )

static

Definition at line 97 of file test_gravity.c.

{
    int i;
    double meanacc = 0, meanforce = 0;
    #pragma omp parallel for reduction(+: meanacc) reduction(+: meanforce)
    for(i = 0; i < PartManager->NumPart; i++)
    {
        int k;
        for(k=0; k<3; k++) {
            if(suppaccns)
                meanacc += fabs(suppaccns[3*i+k]);
            meanforce += fabs(P[i].GravPM[k] + P[i].GravAccel[k]);
        }
    }
    int64_t tot_npart;
    MPI_Allreduce(&PartManager->NumPart, &tot_npart, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    if(suppaccns) {
        MPI_Allreduce(&meanacc, suppmean, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        *suppmean/= (tot_npart*3.);
    }
    MPI_Allreduce(&meanforce, meangrav, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
    *meangrav/= (tot_npart*3.);
}

References GravPM, MPI_INT64, part_manager_type::NumPart, P, and PartManager.

Referenced by check_against_force_direct().

Here is the caller graph for this function:

force_direct()

static void force_direct ( double *  accn )

static

Definition at line 125 of file test_gravity.c.

{
    memset(accn, 0, 3 * sizeof(double) * PartManager->NumPart);
    int xx, yy, zz;
    /* Checked that increasing this has no visible effect on the computed force accuracy*/
    int repeat = 1;
    /* (slowly) compute gravitational force, accounting for periodicity by just inventing extra boxes on either side.*/
    for(xx=-repeat; xx <= repeat; xx++)
        for(yy=-repeat; yy <= repeat; yy++)
            for(zz=-repeat; zz <= repeat; zz++)
            {
                int i;
                double offset[3] = {PartManager->BoxSize * xx, PartManager->BoxSize * yy, PartManager->BoxSize * zz};
                for(i = 0; i < PartManager->NumPart; i++) {
                    int j;
                    for(j = i+1; j < PartManager->NumPart; j++)
                        grav_force(i, j, offset, accn);
                }
            }
}

References part_manager_type::BoxSize, grav_force(), part_manager_type::NumPart, and PartManager.

Referenced by check_against_force_direct().

Here is the call graph for this function:

Here is the caller graph for this function:

grav_force()

static void grav_force ( const int  this, const int  other, const double *  offset, double *  accns  )

static

Definition at line 38 of file test_gravity.c.

{
 
    double r2 = 0;
    int d;
    double dist[3];
    for(d = 0; d < 3; d ++) {
        /* the distance vector points to 'other' */
        dist[d] = offset[d] + P[this].Pos[d] - P[other].Pos[d];
        r2 += dist[d] * dist[d];
    }
 
    const double r = sqrt(r2);
 
    const double h = FORCE_SOFTENING(1, 1);
 
    double fac = 1 / (r2 * r);
    if(r < h) {
        double h_inv = 1.0 / h;
        double h3_inv = h_inv * h_inv * h_inv;
        double u = r * h_inv;
        if(u < 0.5)
            fac = 1. * h3_inv * (10.666666666667 + u * u * (32.0 * u - 38.4));
        else
            fac =
                1. * h3_inv * (21.333333333333 - 48.0 * u +
                        38.4 * u * u - 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
    }
 
    for(d = 0; d < 3; d ++) {
        accns[3*this + d] += - dist[d] * fac * G * P[other].Mass;
        accns[3*other + d] += dist[d] * fac * G * P[this].Mass;
    }
}

References FORCE_SOFTENING(), G, and P.

Referenced by force_direct().

Here is the call graph for this function:

Here is the caller graph for this function:

main()

int main

(

void

)

Definition at line 360 of file test_gravity.c.

               {
    const struct CMUnitTest tests[] = {
        cmocka_unit_test(test_force_flat),
        cmocka_unit_test(test_force_close),
        cmocka_unit_test(test_force_random),
    };
    return cmocka_run_group_tests_mpi(tests, setup_tree, teardown_tree);
}

References setup_tree(), teardown_tree(), test_force_close(), test_force_flat(), and test_force_random().

Here is the call graph for this function:

setup_tree()

static int setup_tree ( void **  state )

static

Definition at line 330 of file test_gravity.c.

                                    {
    walltime_init(&CT);
    /*Set up the important parts of the All structure.*/
    /*Particles should not be outside this*/
    PartManager->BoxSize = 8;
    PartManager->NumPart = 16*16*16;
 
    struct DomainParams dp = {0};
    dp.DomainOverDecompositionFactor = 2;
    dp.DomainUseGlobalSorting = 0;
    dp.TopNodeAllocFactor = 1.;
    dp.SetAsideFactor = 1;
    set_domain_par(dp);
    petapm_module_init(omp_get_max_threads());
    init_forcetree_params(2);
    /*Set up the top-level domain grid*/
    struct forcetree_testdata *data = malloc(sizeof(struct forcetree_testdata));
    data->r = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(data->r, 0);
    *state = (void *) data;
    return 0;
}

References part_manager_type::BoxSize, CT, DomainParams::DomainOverDecompositionFactor, DomainParams::DomainUseGlobalSorting, init_forcetree_params(), part_manager_type::NumPart, PartManager, petapm_module_init(), forcetree_testdata::r, set_domain_par(), DomainParams::SetAsideFactor, DomainParams::TopNodeAllocFactor, and walltime_init().

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 353 of file test_gravity.c.

                                       {
    struct forcetree_testdata * data = (struct forcetree_testdata * ) *state;
    free(data->r);
    free(data);
    return 0;
}

References forcetree_testdata::r.

Referenced by main().

Here is the caller graph for this function:

test_force_close()

static void test_force_close ( void **  state )

static

Definition at line 270 of file test_gravity.c.

                                            {
    /*Set up the particle data*/
    int numpart = PartManager->NumPart;
    int ncbrt = cbrt(numpart);
    double close = 5000;
    P = mymalloc("part", numpart*sizeof(struct particle_data));
    memset(P, 0, 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].Pos[0] = 4. + (i/ncbrt/ncbrt)/close;
        P[i].Pos[1] = 4. + ((i/ncbrt) % ncbrt) /close;
        P[i].Pos[2] = 4. + (i % ncbrt)/close;
    }
    PartManager->NumPart = numpart;
    PartManager->MaxPart = numpart;
    do_force_test(48, 1.5, 0.002, 1);
    myfree(P);
}

References do_force_test(), part_manager_type::MaxPart, myfree, mymalloc, part_manager_type::NumPart, P, and PartManager.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

test_force_flat()

static void test_force_flat ( void **  state )

static

Definition at line 226 of file test_gravity.c.

                                           {
    /*Set up the particle data*/
    int numpart = PartManager->NumPart;
    int ncbrt = cbrt(numpart);
    P = mymalloc("part", numpart*sizeof(struct particle_data));
    memset(P, 0, 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].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);
    }
    PartManager->NumPart = numpart;
    PartManager->MaxPart = numpart;
    do_force_test(48, 1.5, 0.002, 0);
    /* For a homogeneous mass distribution, the force should be zero*/
    double meanerr=0, maxerr=-1;
    #pragma omp parallel for reduction(+: meanerr) reduction(max: maxerr)
    for(i = 0; i < PartManager->NumPart; i++)
    {
        int k;
        for(k=0; k<3; k++) {
            double err = fabs((P[i].GravPM[k] + P[i].GravAccel[k]));
            meanerr += err;
            if(maxerr < err)
                maxerr = err;
        }
    }
    MPI_Allreduce(MPI_IN_PLACE, &meanerr, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(MPI_IN_PLACE, &maxerr, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
    int64_t tot_npart;
    MPI_Allreduce(&PartManager->NumPart, &tot_npart, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
    meanerr/= (tot_npart*3.);
 
    message(0, "Max force %g, mean grav force %g\n", maxerr, meanerr);
    /*Make some statements about the force error*/
    assert_true(maxerr < 0.015);
    assert_true(meanerr < 0.005);
    myfree(P);
}

References part_manager_type::BoxSize, do_force_test(), GravPM, part_manager_type::MaxPart, message(), MPI_INT64, myfree, mymalloc, part_manager_type::NumPart, P, and PartManager.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

test_force_random()

static void test_force_random ( void **  state )

static

Definition at line 316 of file test_gravity.c.

                                             {
    /*Set up the particle data*/
    int numpart = PartManager->NumPart;
    struct forcetree_testdata * data = * (struct forcetree_testdata **) state;
    gsl_rng * r = data->r;
    P = mymalloc("part", numpart*sizeof(struct particle_data));
    memset(P, 0, numpart*sizeof(struct particle_data));
    int i;
    for(i=0; i<2; i++) {
        do_random_test(r, numpart);
    }
    myfree(P);
}

References do_random_test(), myfree, mymalloc, part_manager_type::NumPart, P, PartManager, and forcetree_testdata::r.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

CT

struct ClockTable CT

static

Definition at line 1 of file test_gravity.c.

Referenced by setup_tree().

G

const double G = 43.0071

static

Definition at line 35 of file test_gravity.c.

Referenced by check_omega(), do_force_test(), grav_force(), grav_short_postprocess(), gravpm_init_periodic(), and petapm_init().