runtests.c File Reference

#include <mpi.h>
#include <math.h>
#include "partmanager.h"
#include "forcetree.h"
#include "gravity.h"
#include "petaio.h"
#include "domain.h"
#include "run.h"

Include dependency graph for runtests.c:

Go to the source code of this file.

Functions
char *	GDB_format_particle (int i)
void	register_extra_blocks (struct IOTable *IOTable)
double	copy_and_mean_accn (double(*PairAccn)[3])
void	check_accns (double *meanerr_tot, double *maxerr_tot, double(*PairAccn)[3], double meanacc)
void	run_gravity_test (int RestartSnapNum, Cosmology *CP, const double Asmth, const int Nmesh, const int FastParticleType, const inttime_t Ti_Current, const char *OutputDir, const struct header_data *header)

Function Documentation

check_accns()

void check_accns	(	double *	meanerr_tot,
		double *	maxerr_tot,
		double(*)	PairAccn[3],
		double	meanacc
	)

Definition at line 45 of file runtests.c.

{
    double meanerr=0, maxerr=-1;
    int 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 = 0;
            if(PairAccn[i][k] != 0)
                err = fabs((PairAccn[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 run_gravity_test().

Here is the caller graph for this function:

copy_and_mean_accn()

double copy_and_mean_accn

(

double(*)

PairAccn[3]

)

Definition at line 25 of file runtests.c.

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

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

Referenced by run_gravity_test().

Here is the caller graph for this function:

GDB_format_particle()

char* GDB_format_particle

(

int

i

)

Definition at line 58 of file gdbtools.c.

                                  {
    static char buf[1024];
    char * p = buf;
    int n = 1024;
 
#define add(fmt, ...) \
        snprintf(p, n - 1, fmt, __VA_ARGS__ ); \
        p = buf + strlen(buf); \
        n = 4096 - strlen(buf)
 
    add("P[%d]: ", i);
    add("ID : %ld ", P[i].ID);
    add("Generation: %d ", (int) P[i].Generation);
    add("Mass : %g ", P[i].Mass);
    add("Pos: %g %g %g ", P[i].Pos[0], P[i].Pos[1], P[i].Pos[2]);
    add("Vel: %g %g %g ", P[i].Vel[0], P[i].Vel[1], P[i].Vel[2]);
    add("GravAccel: %g %g %g ", P[i].GravAccel[0], P[i].GravAccel[1], P[i].GravAccel[2]);
    add("GravPM: %g %g %g ", P[i].GravPM[0], P[i].GravPM[1], P[i].GravPM[2]);
    return buf;
}

References add, GravPM, and P.

Referenced by run_gravity_test().

Here is the caller graph for this function:

register_extra_blocks()

void register_extra_blocks

(

struct IOTable *

IOTable

)

Definition at line 16 of file runtests.c.

{
    int ptype;
    for(ptype = 0; ptype < 6; ptype++) {
        IO_REG_WRONLY(GravAccel,       "f4", 3, ptype, IOTable);
        IO_REG_WRONLY(GravPM,       "f4", 3, ptype, IOTable);
    }
}

References GravPM, IO_REG_WRONLY, and ptype.

Referenced by run_gravity_test().

Here is the caller graph for this function:

run_gravity_test()

void run_gravity_test	(	int	RestartSnapNum,
		Cosmology *	CP,
		const double	Asmth,
		const int	Nmesh,
		const int	FastParticleType,
		const inttime_t	Ti_Current,
		const char *	OutputDir,
		const struct header_data *	header
	)

Definition at line 71 of file runtests.c.

{
    DomainDecomp ddecomp[1] = {0};
    domain_decompose_full(ddecomp);
 
    struct IOTable IOTable = {0};
    /* NO metals written*/
    register_io_blocks(&IOTable, 0, 0);
    register_extra_blocks(&IOTable);
 
    double (* PairAccn)[3] = (double (*) [3]) mymalloc2("PairAccns", 3*sizeof(double) * PartManager->NumPart);
 
    PetaPM pm[1] = {0};
    gravpm_init_periodic(pm, PartManager->BoxSize, Asmth, Nmesh, CP->GravInternal);
 
    int NTask;
    MPI_Comm_size(MPI_COMM_WORLD, &NTask);
    ActiveParticles Act = {0};
    DriftKickTimes times = init_driftkicktime(Ti_Current);
    rebuild_activelist(&Act, &times, 0, header->TimeSnapshot);
 
    ForceTree Tree = {0};
    force_tree_rebuild(&Tree, ddecomp, 1, 1, OutputDir);
    gravpm_force(pm, &Tree, CP, header->TimeSnapshot, header->UnitLength_in_cm, OutputDir, header->TimeIC, FastParticleType);
    force_tree_rebuild(&Tree, ddecomp, 1, 1, OutputDir);
 
    struct gravshort_tree_params origtreeacc = get_gravshort_treepar();
    struct gravshort_tree_params treeacc = origtreeacc;
    const double rho0 = CP->Omega0 * 3 * CP->Hubble * CP->Hubble / (8 * M_PI * CP->GravInternal);
    grav_short_pair(&Act, pm, &Tree, treeacc.Rcut, rho0, 0, FastParticleType);
 
    double meanacc = copy_and_mean_accn(PairAccn);
    message(0, "GravShort Pairs %s\n", GDB_format_particle(0));
    char * fname = fastpm_strdup_printf("%s/PART-pairs-%03d", OutputDir, RestartSnapNum);
 
    petaio_save_snapshot(fname, &IOTable, 0, header->TimeSnapshot, CP);
 
    treeacc.ErrTolForceAcc = 0;
    set_gravshort_treepar(treeacc);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
 
    /* This checks fully opened tree force against pair force*/
    double meanerr, maxerr;
    check_accns(&meanerr,&maxerr,PairAccn, meanacc);
    message(0, "Force error, open tree vs pairwise. max : %g mean: %g forcetol: %g\n", maxerr, meanerr, treeacc.ErrTolForceAcc);
 
    if(maxerr > 0.1)
        endrun(2, "Fully open tree force does not agree with pairwise calculation! maxerr %g > 0.1!\n", maxerr);
 
    message(0, "GravShort Tree %s\n", GDB_format_particle(0));
    fname = fastpm_strdup_printf("%s/PART-tree-open-%03d", OutputDir, RestartSnapNum);
    petaio_save_snapshot(fname, &IOTable, 0, header->TimeSnapshot, CP);
 
    /* This checks tree force against tree force with zero error (which always opens).*/
    copy_and_mean_accn(PairAccn);
 
    treeacc = origtreeacc;
    set_gravshort_treepar(treeacc);
    /* Code automatically sets the UseTreeBH parameter.*/
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
 
    fname = fastpm_strdup_printf("%s/PART-tree-%03d", OutputDir, RestartSnapNum);
    petaio_save_snapshot(fname, &IOTable, 0, header->TimeSnapshot, CP);
 
    check_accns(&meanerr,&maxerr,PairAccn, meanacc);
    message(0, "Force error, open tree vs tree.: %g mean: %g forcetol: %g\n", maxerr, meanerr, treeacc.ErrTolForceAcc);
 
    if(meanerr > treeacc.ErrTolForceAcc* 1.2)
        endrun(2, "Average force error is underestimated: %g > 1.2 * %g!\n", meanerr, treeacc.ErrTolForceAcc);
 
    copy_and_mean_accn(PairAccn);
    /* This checks the tree against a larger Rcut.*/
    treeacc.Rcut = 9.5;
    set_gravshort_treepar(treeacc);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
    fname = fastpm_strdup_printf("%s/PART-tree-rcut-%03d", OutputDir, RestartSnapNum);
    petaio_save_snapshot(fname, &IOTable, 0, header->TimeSnapshot, CP);
 
    check_accns(&meanerr,&maxerr,PairAccn, meanacc);
    message(0, "Force error, tree vs rcut.: %g mean: %g Rcut = %g\n", maxerr, meanerr, treeacc.Rcut);
 
    if(maxerr > 0.2)
        endrun(2, "Rcut decreased below desired value, error too large %g\n", maxerr);
 
    /* This checks the tree against a box with a smaller Nmesh.*/
    treeacc = origtreeacc;
    force_tree_free(&Tree);
 
    petapm_destroy(pm);
 
    gravpm_init_periodic(pm, PartManager->BoxSize, Asmth, Nmesh/2., CP->GravInternal);
    force_tree_rebuild(&Tree, ddecomp, 1, 1, OutputDir);
    gravpm_force(pm, &Tree, CP, header->TimeSnapshot, header->UnitLength_in_cm, OutputDir, header->TimeIC, FastParticleType);
    force_tree_rebuild(&Tree, ddecomp, 1, 1, OutputDir);
    set_gravshort_treepar(treeacc);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
    grav_short_tree(&Act, pm, &Tree, rho0, 0, FastParticleType);
    fname = fastpm_strdup_printf("%s/PART-tree-nmesh2-%03d", OutputDir, RestartSnapNum);
    petaio_save_snapshot(fname, &IOTable, 0, header->TimeSnapshot, CP);
 
    check_accns(&meanerr, &maxerr, PairAccn, meanacc);
    message(0, "Force error, nmesh %d vs %d: %g mean: %g \n", Nmesh, Nmesh/2, maxerr, meanerr);
 
    if(maxerr > 0.5 || meanerr > 0.05)
        endrun(2, "Nmesh sensitivity worse, something may be wrong\n");
 
    force_tree_free(&Tree);
    petapm_destroy(pm);
 
    myfree(PairAccn);
 
    destroy_io_blocks(&IOTable);
    domain_free(ddecomp);
}

References part_manager_type::BoxSize, check_accns(), copy_and_mean_accn(), CP, destroy_io_blocks(), domain_decompose_full(), domain_free(), endrun(), gravshort_tree_params::ErrTolForceAcc, fastpm_strdup_printf(), force_tree_free(), force_tree_rebuild(), GDB_format_particle(), get_gravshort_treepar(), grav_short_pair(), grav_short_tree(), Cosmology::GravInternal, gravpm_force(), gravpm_init_periodic(), Cosmology::Hubble, init_driftkicktime(), message(), myfree, mymalloc2, NTask, part_manager_type::NumPart, Cosmology::Omega0, PartManager, petaio_save_snapshot(), petapm_destroy(), gravshort_tree_params::Rcut, rebuild_activelist(), register_extra_blocks(), register_io_blocks(), set_gravshort_treepar(), header_data::TimeIC, header_data::TimeSnapshot, and header_data::UnitLength_in_cm.

Referenced by runtests().

Here is the call graph for this function:

Here is the caller graph for this function: