gravshort-tree.c File Reference

#include <mpi.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <sys/ipc.h>
#include <sys/sem.h>
#include "utils.h"
#include "forcetree.h"
#include "treewalk.h"
#include "timestep.h"
#include "gravshort.h"
#include "walltime.h"

Include dependency graph for gravshort-tree.c:

Go to the source code of this file.

Functions
double	FORCE_SOFTENING (int i, int type)
void	gravshort_set_softenings (double MeanSeparation)
void	set_gravshort_treepar (struct gravshort_tree_params tree_params)
struct gravshort_tree_params	get_gravshort_treepar (void)
void	set_gravshort_tree_params (ParameterSet *ps)
int	force_treeev_shortrange (TreeWalkQueryGravShort *input, TreeWalkResultGravShort *output, LocalTreeWalk *lv)
void	grav_short_tree (const ActiveParticles *act, PetaPM *pm, ForceTree *tree, double rho0, int NeutrinoTracer, int FastParticleType)
static void	apply_accn_to_output (TreeWalkResultGravShort *output, const double dx[3], const double r2, const double h, const double mass, const double cellsize)
static int	shall_we_discard_node (const double len, const double r2, const double center[3], const double inpos[3], const double BoxSize, const double rcut, const double rcut2)
static int	shall_we_open_node (const double len, const double mass, const double r2, const double center[3], const double inpos[3], const double BoxSize, const double aold, const int TreeUseBH, const double BHOpeningAngle2)

Variables
static struct gravshort_tree_params	TreeParams
double	GravitySoftening

Function Documentation

apply_accn_to_output()

static void apply_accn_to_output ( TreeWalkResultGravShort *  output, const double  dx[3], const double  r2, const double  h, const double  mass, const double  cellsize  )

static

Definition at line 176 of file gravshort-tree.c.

{
    const double r = sqrt(r2);
 
    double fac = mass / (r2 * r);
    double facpot = -mass / r;
 
    if(r2 < h*h)
    {
        double wp;
        const double h3_inv = 1.0 / h / h / h;
        const double u = r / h;
        if(u < 0.5) {
            fac = mass * h3_inv * (10.666666666667 + u * u * (32.0 * u - 38.4));
            wp = -2.8 + u * u * (5.333333333333 + u * u * (6.4 * u - 9.6));
        }
        else {
            fac =
                mass * h3_inv * (21.333333333333 - 48.0 * u +
                        38.4 * u * u - 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
            wp =
                -3.2 + 0.066666666667 / u + u * u * (10.666666666667 +
                        u * (-16.0 + u * (9.6 - 2.133333333333 * u)));
        }
        facpot = mass / h * wp;
    }
 
    if(0 == grav_apply_short_range_window(r, &fac, &facpot, cellsize)) {
        int i;
        for(i = 0; i < 3; i++)
            output->Acc[i] += dx[i] * fac;
        output->Potential += facpot;
    }
}

References TreeWalkResultGravShort::Acc, GravShortPriv::cellsize, grav_apply_short_range_window(), and TreeWalkResultGravShort::Potential.

Referenced by force_treeev_shortrange().

Here is the call graph for this function:

Here is the caller graph for this function:

FORCE_SOFTENING()

double FORCE_SOFTENING	(	int	i,
		int	type
	)

Definition at line 36 of file gravshort-tree.c.

{
    if (TreeParams.AdaptiveSoftening == 1 && type == 0) {
        return P[i].Hsml;
    }
    /* Force is Newtonian beyond this.*/
    return 2.8 * GravitySoftening;
}

References gravshort_tree_params::AdaptiveSoftening, GravitySoftening, P, and TreeParams.

Referenced by blackhole_accretion_ngbiter(), density(), force_treeev_shortrange(), get_timestep_dloga(), grav_force(), grav_short_copy(), grav_short_pair_ngbiter(), and grav_short_postprocess().

Here is the caller graph for this function:

force_treeev_shortrange()

int force_treeev_shortrange	(	TreeWalkQueryGravShort *	input,
		TreeWalkResultGravShort *	output,
		LocalTreeWalk *	lv
	)

In the TreePM algorithm, the tree is walked only locally around the target coordinate. Tree nodes that fall outside a box of half side-length Rcut= RCUT*ASMTH*MeshSize can be discarded. The short-range potential is modified by a complementary error function, multiplied with the Newtonian form. The resulting short-range suppression compared to the Newtonian force is tabulated, because looking up from this table is faster than recomputing the corresponding factor, despite the memory-access penalty (which reduces cache performance) incurred by the table.

Definition at line 267 of file gravshort-tree.c.

{
    const ForceTree * tree = lv->tw->tree;
    const double BoxSize = tree->BoxSize;
 
    /*Tree-opening constants*/
    const double cellsize = GRAV_GET_PRIV(lv->tw)->cellsize;
    const double rcut = GRAV_GET_PRIV(lv->tw)->Rcut;
    const double rcut2 = rcut * rcut;
    const double aold = GRAV_GET_PRIV(lv->tw)->ErrTolForceAcc * input->OldAcc;
    const int TreeUseBH = GRAV_GET_PRIV(lv->tw)->TreeUseBH;
    const double BHOpeningAngle2 = GRAV_GET_PRIV(lv->tw)->BHOpeningAngle * GRAV_GET_PRIV(lv->tw)->BHOpeningAngle;
    const int NeutrinoTracer = GRAV_GET_PRIV(lv->tw)->NeutrinoTracer;
    const int FastParticleType = GRAV_GET_PRIV(lv->tw)->FastParticleType;
 
    /*Input particle data*/
    const double * inpos = input->base.Pos;
 
    /*Start the tree walk*/
    int listindex;
 
    /* Primary treewalk only ever has one nodelist entry*/
    for(listindex = 0; listindex < NODELISTLENGTH && (lv->mode == 1 || listindex < 1); listindex++)
    {
        int numcand = 0;
        /* Use the next node in the node list if we are doing a secondary walk.
         * For a primary walk the node list only ever contains one node. */
        int no = input->base.NodeList[listindex];
        int startno = no;
        if(no < 0)
            break;
 
        while(no >= 0)
        {
            /* The tree always walks internal nodes*/
            struct NODE *nop = &tree->Nodes[no];
 
            if(lv->mode == 1)
            {
                if(nop->f.TopLevel && no != startno)    /* we reached a top-level node again, which means that we are done with the branch */
                {
                    no = -1;
                    continue;
                }
            }
 
            int i;
            double dx[3];
            for(i = 0; i < 3; i++)
                dx[i] = NEAREST(nop->mom.cofm[i] - inpos[i], BoxSize);
            const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
 
            /* Discard this node, move to sibling*/
            if(shall_we_discard_node(nop->len, r2, nop->center, inpos, BoxSize, rcut, rcut2))
            {
                no = nop->sibling;
                /* Don't add this node*/
                continue;
            }
 
            /* This node accelerates the particle directly, and is not opened.*/
            int open_node = shall_we_open_node(nop->len, nop->mom.mass, r2, nop->center, inpos, BoxSize, aold, TreeUseBH, BHOpeningAngle2);
            if(TreeParams.AdaptiveSoftening == 1 && (input->Soft < nop->mom.hmax))
            {
                /* Always open the node if it has a larger softening than the particle,
                 * and the particle is inside its softening radius.
                 * This condition only ever applies for adaptive softenings. It may or may not make sense. */
                if(r2 < nop->mom.hmax * nop->mom.hmax)
                    open_node = 1;
            }
 
            if(!open_node)
            {
                double h = input->Soft;
                if(TreeParams.AdaptiveSoftening)
                    h = DMAX(input->Soft, nop->mom.hmax);
                /* ok, node can be used */
                no = nop->sibling;
                /* Compute the acceleration and apply it to the output structure*/
                apply_accn_to_output(output, dx, r2, h, nop->mom.mass, cellsize);
                continue;
            }
 
            /* Now we have a cell that needs to be opened.
             * If it contains particles we can add them directly here */
            if(nop->f.ChildType == PARTICLE_NODE_TYPE)
            {
                /* Loop over child particles*/
                for(i = 0; i < nop->s.noccupied; i++) {
                    int pp = nop->s.suns[i];
                    lv->ngblist[numcand++] = pp;
                }
                no = nop->sibling;
            }
            else if (nop->f.ChildType == PSEUDO_NODE_TYPE)
            {
                if(lv->mode == 0)
                {
                    if(-1 == treewalk_export_particle(lv, nop->s.suns[0]))
                        return -1;
                }
 
                /* Move to the sibling (likely also a pseudo node)*/
                no = nop->sibling;
            }
            else if(nop->f.ChildType == NODE_NODE_TYPE)
            {
                /* This node contains other nodes and we need to open it.*/
                no = nop->s.suns[0];
            }
        }
        int i;
        for(i = 0; i < numcand; i++)
        {
            int pp = lv->ngblist[i];
            /* Fast particle neutrinos don't cause short-range acceleration before activation.*/
            if(NeutrinoTracer && P[pp].Type == FastParticleType)
                continue;
 
            double dx[3];
            int j;
            for(j = 0; j < 3; j++)
                dx[j] = NEAREST(P[pp].Pos[j] - inpos[j], BoxSize);
            const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
 
            /* This is always the Newtonian softening,
             * match the default from FORCE_SOFTENING. */
            double h = 2.8 * GravitySoftening;
            if(TreeParams.AdaptiveSoftening == 1) {
                h = DMAX(input->Soft, FORCE_SOFTENING(pp, P[pp].Type));
            }
            /* Compute the acceleration and apply it to the output structure*/
            apply_accn_to_output(output, dx, r2, h, P[pp].Mass, cellsize);
        }
        lv->Ninteractions += numcand;
    }
 
    if(lv->mode == 1) {
        lv->Nnodesinlist += listindex;
        lv->Nlist += 1;
    }
    return 1;
}

References gravshort_tree_params::AdaptiveSoftening, apply_accn_to_output(), TreeWalkQueryGravShort::base, ForceTree::BoxSize, GravShortPriv::cellsize, NODE::center, NODE::ChildType, NODE::cofm, DMAX, NODE::f, GravShortPriv::FastParticleType, FORCE_SOFTENING(), GRAV_GET_PRIV, GravitySoftening, NODE::hmax, NODE::len, NODE::mass, LocalTreeWalk::mode, NODE::mom, NEAREST, GravShortPriv::NeutrinoTracer, LocalTreeWalk::ngblist, LocalTreeWalk::Ninteractions, LocalTreeWalk::Nlist, LocalTreeWalk::Nnodesinlist, NodeChild::noccupied, NODE_NODE_TYPE, TreeWalkQueryBase::NodeList, NODELISTLENGTH, ForceTree::Nodes, TreeWalkQueryGravShort::OldAcc, P, PARTICLE_NODE_TYPE, TreeWalkQueryBase::Pos, PSEUDO_NODE_TYPE, NODE::s, shall_we_discard_node(), shall_we_open_node(), NODE::sibling, TreeWalkQueryGravShort::Soft, NodeChild::suns, NODE::TopLevel, TreeWalk::tree, TreeParams, GravShortPriv::TreeUseBH, treewalk_export_particle(), and LocalTreeWalk::tw.

Referenced by grav_short_tree().

Here is the call graph for this function:

Here is the caller graph for this function:

get_gravshort_treepar()

struct gravshort_tree_params get_gravshort_treepar

(

void

)

Definition at line 55 of file gravshort-tree.c.

{
    return TreeParams;
}

References TreeParams.

Referenced by run(), and run_gravity_test().

Here is the caller graph for this function:

grav_short_tree()

void grav_short_tree	(	const ActiveParticles *	act,
		PetaPM *	pm,
		ForceTree *	tree,
		double	rho0,
		int	NeutrinoTracer,
		int	FastParticleType
	)

This function computes the gravitational forces for all active particles. If needed, a new tree is constructed, otherwise the dynamically updated tree is used. Particles are only exported to other processors when really needed, thereby allowing a good use of the communication buffer. NeutrinoTracer = All.HybridNeutrinosOn && (atime <= All.HybridNuPartTime); rho0 = CP.Omega0 * 3 * CP.Hubble * CP.Hubble / (8 * M_PI * G)

Definition at line 105 of file gravshort-tree.c.

{
    double timeall = 0;
    double timetree, timewait, timecomm;
 
    TreeWalk tw[1] = {{0}};
    struct GravShortPriv priv;
    priv.cellsize = tree->BoxSize / pm->Nmesh;
    priv.Rcut = TreeParams.Rcut * pm->Asmth * priv.cellsize;;
    priv.ErrTolForceAcc = TreeParams.ErrTolForceAcc;
    priv.TreeUseBH = TreeParams.TreeUseBH;
    priv.BHOpeningAngle = TreeParams.BHOpeningAngle;
    priv.FastParticleType = FastParticleType;
    priv.NeutrinoTracer = NeutrinoTracer;
    priv.G = pm->G;
    priv.cbrtrho0 = pow(rho0, 1.0 / 3);
 
    if(!tree->moments_computed_flag)
        endrun(2, "Gravtree called before tree moments computed!\n");
 
    tw->ev_label = "GRAVTREE";
    tw->visit = (TreeWalkVisitFunction) force_treeev_shortrange;
    /* gravity applies to all particles. Including Tracer particles to enhance numerical stability. */
    tw->haswork = NULL;
    tw->reduce = (TreeWalkReduceResultFunction) grav_short_reduce;
    tw->postprocess = (TreeWalkProcessFunction) grav_short_postprocess;
 
    tw->query_type_elsize = sizeof(TreeWalkQueryGravShort);
    tw->result_type_elsize = sizeof(TreeWalkResultGravShort);
    tw->fill = (TreeWalkFillQueryFunction) grav_short_copy;
    tw->tree = tree;
    tw->priv = &priv;
 
    walltime_measure("/Misc");
 
    /* allocate buffers to arrange communication */
    MPIU_Barrier(MPI_COMM_WORLD);
    message(0, "Begin tree force.  (presently allocated=%g MB)\n", mymalloc_usedbytes() / (1024.0 * 1024.0));
 
    walltime_measure("/Misc");
 
    treewalk_run(tw, act->ActiveParticle, act->NumActiveParticle);
 
    /* Now the force computation is finished */
    /*  gather some diagnostic information */
 
    timetree = tw->timecomp1 + tw->timecomp2 + tw->timecomp3;
    timewait = tw->timewait1 + tw->timewait2;
    timecomm= tw->timecommsumm1 + tw->timecommsumm2;
 
    walltime_add("/Tree/Walk1", tw->timecomp1);
    walltime_add("/Tree/Walk2", tw->timecomp2);
    walltime_add("/Tree/PostProcess", tw->timecomp3);
    walltime_add("/Tree/Send", tw->timecommsumm1);
    walltime_add("/Tree/Recv", tw->timecommsumm2);
    walltime_add("/Tree/Wait1", tw->timewait1);
    walltime_add("/Tree/Wait2", tw->timewait2);
 
    timeall = walltime_measure(WALLTIME_IGNORE);
 
    walltime_add("/Tree/Misc", timeall - (timetree + timewait + timecomm));
 
    /* TreeUseBH > 1 means use the BH criterion on the initial timestep only,
     * avoiding the fully open O(N^2) case.*/
    if(TreeParams.TreeUseBH > 1)
        TreeParams.TreeUseBH = 0;
}

References ActiveParticles::ActiveParticle, PetaPM::Asmth, gravshort_tree_params::BHOpeningAngle, GravShortPriv::BHOpeningAngle, ForceTree::BoxSize, GravShortPriv::cbrtrho0, GravShortPriv::cellsize, endrun(), gravshort_tree_params::ErrTolForceAcc, GravShortPriv::ErrTolForceAcc, TreeWalk::ev_label, GravShortPriv::FastParticleType, TreeWalk::fill, force_treeev_shortrange(), GravShortPriv::G, PetaPM::G, grav_short_copy(), grav_short_postprocess(), grav_short_reduce(), TreeWalk::haswork, message(), ForceTree::moments_computed_flag, MPIU_Barrier, mymalloc_usedbytes, GravShortPriv::NeutrinoTracer, PetaPM::Nmesh, ActiveParticles::NumActiveParticle, TreeWalk::postprocess, TreeWalk::priv, TreeWalk::query_type_elsize, gravshort_tree_params::Rcut, GravShortPriv::Rcut, TreeWalk::reduce, TreeWalk::result_type_elsize, TreeWalk::timecommsumm1, TreeWalk::timecommsumm2, TreeWalk::timecomp1, TreeWalk::timecomp2, TreeWalk::timecomp3, TreeWalk::timewait1, TreeWalk::timewait2, TreeWalk::tree, TreeParams, gravshort_tree_params::TreeUseBH, GravShortPriv::TreeUseBH, treewalk_run(), TreeWalk::visit, walltime_add, WALLTIME_IGNORE, and walltime_measure.

Referenced by do_force_test(), run(), and run_gravity_test().

Here is the call graph for this function:

Here is the caller graph for this function:

gravshort_set_softenings()

void gravshort_set_softenings

(

double

MeanSeparation

)

Sets the (comoving) softening length, converting from units of the mean separation to comoving internal units.

Definition at line 47 of file gravshort-tree.c.

{
    GravitySoftening = TreeParams.FractionalGravitySoftening * MeanSeparation;
    /* 0: Gas is collisional */
    message(0, "GravitySoftening = %g\n", GravitySoftening);
}

References gravshort_tree_params::FractionalGravitySoftening, GravitySoftening, message(), and TreeParams.

Referenced by do_force_test(), init(), and setup_density().

Here is the call graph for this function:

Here is the caller graph for this function:

set_gravshort_tree_params()

void set_gravshort_tree_params

(

ParameterSet *

ps

)

Definition at line 67 of file gravshort-tree.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        TreeParams.BHOpeningAngle = param_get_double(ps, "BHOpeningAngle");
        TreeParams.ErrTolForceAcc = param_get_double(ps, "ErrTolForceAcc");
        TreeParams.BHOpeningAngle = param_get_double(ps, "BHOpeningAngle");
        TreeParams.TreeUseBH= param_get_int(ps, "TreeUseBH");
        TreeParams.Rcut = param_get_double(ps, "TreeRcut");
        TreeParams.FractionalGravitySoftening = param_get_double(ps, "GravitySoftening");
        TreeParams.AdaptiveSoftening = !param_get_int(ps, "GravitySofteningGas");
 
 
    }
    MPI_Bcast(&TreeParams, sizeof(struct gravshort_tree_params), MPI_BYTE, 0, MPI_COMM_WORLD);
}

References gravshort_tree_params::AdaptiveSoftening, gravshort_tree_params::BHOpeningAngle, gravshort_tree_params::ErrTolForceAcc, gravshort_tree_params::FractionalGravitySoftening, param_get_double(), param_get_int(), gravshort_tree_params::Rcut, ThisTask, TreeParams, and gravshort_tree_params::TreeUseBH.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

set_gravshort_treepar()

void set_gravshort_treepar

(

struct gravshort_tree_params

tree_params

)

Definition at line 55 of file gravshort-tree.c.

{
    TreeParams = tree_params;
}

Referenced by do_force_test(), run_gravity_test(), and setup_density().

Here is the caller graph for this function:

shall_we_discard_node()

static int shall_we_discard_node ( const double  len, const double  r2, const double  center[3], const double  inpos[3], const double  BoxSize, const double  rcut, const double  rcut2  )

static

Definition at line 215 of file gravshort-tree.c.

{
    /* This checks the distance from the node center of mass
     * is greater than the cutoff. */
    if(r2 > rcut2)
    {
        /* check whether we can stop walking along this branch */
        const double eff_dist = rcut + 0.5 * len;
        int i;
        /*This checks whether we are also outside this region of the oct-tree*/
        /* As long as one dimension is outside, we are fine*/
        for(i=0; i < 3; i++)
            if(fabs(NEAREST(center[i] - inpos[i], BoxSize)) > eff_dist)
                return 1;
    }
    return 0;
}

References NEAREST.

Referenced by force_treeev_shortrange().

Here is the caller graph for this function:

shall_we_open_node()

static int shall_we_open_node ( const double  len, const double  mass, const double  r2, const double  center[3], const double  inpos[3], const double  BoxSize, const double  aold, const int  TreeUseBH, const double  BHOpeningAngle2  )

static

Definition at line 237 of file gravshort-tree.c.

{
    /* Check the relative acceleration opening condition*/
    if((TreeUseBH == 0) && (mass * len * len > r2 * r2 * aold))
         return 1;
     /*Check Barnes-Hut opening angle*/
    if((TreeUseBH > 0) && (len * len > r2 * BHOpeningAngle2))
         return 1;
 
    const double inside = 0.6 * len;
    /* Open the cell if we are inside it, even if the opening criterion is not satisfied.*/
    if(fabs(NEAREST(center[0] - inpos[0], BoxSize)) < inside &&
        fabs(NEAREST(center[1] - inpos[1], BoxSize)) < inside &&
        fabs(NEAREST(center[2] - inpos[2], BoxSize)) < inside)
        return 1;
 
    /* ok, node can be used */
    return 0;
}

References NEAREST, and GravShortPriv::TreeUseBH.

Referenced by force_treeev_shortrange().

Here is the caller graph for this function:

Variable Documentation

GravitySoftening

double GravitySoftening

Definition at line 31 of file gravshort-tree.c.

Referenced by FORCE_SOFTENING(), force_treeev_shortrange(), and gravshort_set_softenings().

TreeParams

struct gravshort_tree_params TreeParams

static

Definition at line 1 of file gravshort-tree.c.

Referenced by FORCE_SOFTENING(), force_treeev_shortrange(), get_gravshort_treepar(), grav_short_tree(), gravshort_set_softenings(), and set_gravshort_tree_params().