forcetree.c File Reference

gravitational tree More...

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <omp.h>
#include "slotsmanager.h"
#include "partmanager.h"
#include "domain.h"
#include "forcetree.h"
#include "checkpoint.h"
#include "walltime.h"
#include "utils/endrun.h"

Include dependency graph for forcetree.c:

Go to the source code of this file.

Classes
struct	forcetree_params
struct	NodeCache

Macros
#define	NODECACHE_SIZE   (8*12)

Functions
void	init_forcetree_params (const int FastParticleType)
static ForceTree	force_tree_build (int npart, int mask, DomainDecomp *ddecomp, const int HybridNuGrav, const int DoMoments, const char *EmergencyOutputDir)
static void	force_treeupdate_pseudos (int no, const ForceTree *tree)
static void	force_create_node_for_topnode (int no, int topnode, struct NODE *Nodes, const DomainDecomp *ddecomp, int bits, int x, int y, int z, int *nextfree, const int lastnode)
static void	force_exchange_pseudodata (ForceTree *tree, const DomainDecomp *ddecomp)
static void	force_insert_pseudo_particles (const ForceTree *tree, const DomainDecomp *ddecomp)
static void	add_particle_moment_to_node (struct NODE *pnode, int i)
static int	force_tree_eh_slots_fork (EIBase *event, void *userdata)
int	force_tree_allocated (const ForceTree *tree)
void	force_tree_rebuild (ForceTree *tree, DomainDecomp *ddecomp, const int HybridNuGrav, const int DoMoments, const char *EmergencyOutputDir)
void	force_tree_rebuild_mask (ForceTree *tree, DomainDecomp *ddecomp, int mask, const int HybridNuGrav, const char *EmergencyOutputDir)
int	get_subnode (const struct NODE *node, const int p_i)
static int	inside_node (const struct NODE *node, const int p_i)
static void	init_internal_node (struct NODE *nfreep, struct NODE *parent, int subnode)
int	get_freenode (int *nnext, struct NodeCache *nc)
static int	modify_internal_node (int parent, int subnode, int p_toplace, const ForceTree tb, const int HybridNuGrav)
static int	create_new_node_layer (int firstparent, int p_toplace, const int HybridNuGrav, const ForceTree tb, int *nnext, struct NodeCache *nc)
int	add_particle_to_tree (int i, int cur_start, const ForceTree tb, const int HybridNuGrav, struct NodeCache *nc, int *nnext)
int	merge_partial_force_trees (int left, int right, struct NodeCache *nc, int *nnext, const struct ForceTree tb, int HybridNuGrav)
int	force_tree_create_nodes (const ForceTree tb, const int npart, int mask, DomainDecomp *ddecomp, const int HybridNuGrav)
int	force_get_father (int no, const ForceTree *tree)
static int	force_get_sibling (const int sib, const int j, const int *suns)
static void	force_update_particle_node (int no, const ForceTree *tree)
static int	force_update_node_recursive (int no, int sib, int level, const ForceTree *tree)
void	force_update_node_parallel (const ForceTree *tree, const DomainDecomp *ddecomp)
void	force_update_hmax (int *activeset, int size, ForceTree *tree, DomainDecomp *ddecomp)
ForceTree	force_treeallocate (int64_t maxnodes, int64_t maxpart, DomainDecomp *ddecomp)
void	force_tree_free (ForceTree *tree)

Variables
static struct forcetree_params	ForceTreeParams

Detailed Description

gravitational tree

This file contains the computation of the gravitational force by means of a tree. The type of tree implemented is a geometrical oct-tree, starting from a cube encompassing all particles. This cube is automatically found in the ddecomp decomposition, which also splits up the global "top-level" tree along node boundaries, moving the particles of different parts of the tree to separate processors.

Local naming convention: once built it is ForceTree * tree, passed by reference. While the nodes are still being added it is ForceTree tb, passed by value.

Definition in file forcetree.c.

Macro Definition Documentation

NODECACHE_SIZE

#define NODECACHE_SIZE   (8*12)

Definition at line 336 of file forcetree.c.

Function Documentation

add_particle_moment_to_node()

static void add_particle_moment_to_node ( struct NODE *  pnode, int  i  )

static

Definition at line 914 of file forcetree.c.

{
    int k;
    pnode->mom.mass += (P[i].Mass);
    for(k=0; k<3; k++)
        pnode->mom.cofm[k] += (P[i].Mass * P[i].Pos[k]);
 
    if(P[i].Type == 0)
    {
        int j;
        /* Maximal distance any of the member particles peek out from the side of the node.
         * May be at most hmax, as |Pos - Center| < len.*/
        for(j = 0; j < 3; j++) {
            pnode->mom.hmax = DMAX(pnode->mom.hmax, fabs(P[i].Pos[j] - pnode->center[j]) + P[i].Hsml - pnode->len);
        }
    }
}

References NODE::center, NODE::cofm, DMAX, NODE::hmax, NODE::len, NODE::mass, NODE::mom, and P.

Referenced by modify_internal_node().

Here is the caller graph for this function:

add_particle_to_tree()

int add_particle_to_tree	(	int	i,
		int	cur_start,
		const ForceTree	tb,
		const int	HybridNuGrav,
		struct NodeCache *	nc,
		int *	nnext
	)

Definition at line 487 of file forcetree.c.

{
    int child, nocc;
    int cur = cur_start;
    /*Walk the main tree until we get something that isn't an internal node.*/
    do
    {
        /*No lock needed: if we have an internal node here it will be stable*/
        nocc = tb.Nodes[cur].s.noccupied;
 
        /* This node still has space for a particle (or needs conversion)*/
        if(nocc < (1 << 16))
            break;
 
        /* This node has child subnodes: find them.*/
        int subnode = get_subnode(&tb.Nodes[cur], i);
        /*No lock needed: if we have an internal node here it will be stable*/
        child = tb.Nodes[cur].s.suns[subnode];
 
        if(child > tb.lastnode || child < tb.firstnode)
            endrun(1,"Corruption in tree build: N[%d].[%d] = %d > lastnode (%d)\n",cur, subnode, child, tb.lastnode);
        cur = child;
    }
    while(child >= tb.firstnode);
 
    /* We have a guaranteed spot.*/
    nocc = tb.Nodes[cur].s.noccupied;
    tb.Nodes[cur].s.noccupied++;
 
    /* Now we have something that isn't an internal node. We can place the particle! */
    if(nocc < NMAXCHILD)
        modify_internal_node(cur, nocc, i, tb, HybridNuGrav);
    /* In this case we need to create a new layer of nodes beneath this one*/
    else if(nocc < 1<<16) {
        if(create_new_node_layer(cur, i, HybridNuGrav, tb, nnext, nc))
            return -1;
    } else
        endrun(2, "Tried to convert already converted node %d with nocc = %d\n", cur, nocc);
    return cur;
}

References create_new_node_layer(), endrun(), ForceTree::firstnode, get_subnode(), ForceTree::lastnode, modify_internal_node(), NMAXCHILD, NodeChild::noccupied, ForceTree::Nodes, NODE::s, and NodeChild::suns.

Referenced by force_tree_create_nodes(), and merge_partial_force_trees().

Here is the call graph for this function:

Here is the caller graph for this function:

create_new_node_layer()

static int create_new_node_layer ( int  firstparent, int  p_toplace, const int  HybridNuGrav, const ForceTree  tb, int *  nnext, struct NodeCache *  nc  )

static

Definition at line 376 of file forcetree.c.

{
    /* This is so we can defer changing
     * the type of the existing node until the end.*/
    int parent = firstparent;
 
    do {
        int i;
        struct NODE *nprnt = &tb.Nodes[parent];
 
        /* Braces to scope oldsuns and newsuns*/
        {
        int newsuns[NMAXCHILD];
 
        int * oldsuns = nprnt->s.suns;
 
        /*We have two particles here, so create a new child node to store them both.*/
        /* if we are here the node must be large enough, thus contain exactly one child. */
        /* The parent is already a leaf, need to split */
        /* Get memory for 8 extra nodes from our cache.*/
        newsuns[0] = get_freenode(nnext, nc);
        /*If we already have too many nodes, exit loop.*/
        if(nc->nnext_thread >= tb.lastnode) {
            /* This means that we have > NMAXCHILD particles in the same place,
            * which usually indicates a bug in the particle evolution. Print some helpful debug information.*/
            message(1, "Failed placing %d at %g %g %g, type %d, ID %ld. Others were %d (%g %g %g, t %d ID %ld) and %d (%g %g %g, t %d ID %ld). next %d last %d\n",
                p_toplace, P[p_toplace].Pos[0], P[p_toplace].Pos[1], P[p_toplace].Pos[2], P[p_toplace].Type, P[p_toplace].ID,
                oldsuns[0], P[oldsuns[0]].Pos[0], P[oldsuns[0]].Pos[1], P[oldsuns[0]].Pos[2], P[oldsuns[0]].Type, P[oldsuns[0]].ID,
                oldsuns[1], P[oldsuns[1]].Pos[0], P[oldsuns[1]].Pos[1], P[oldsuns[1]].Pos[2], P[oldsuns[1]].Type, P[oldsuns[1]].ID
            );
            nc->nnext_thread = tb.lastnode + 10 * NODECACHE_SIZE;
            /* If this is not the first layer created,
                * we need to mark the overall parent as a node node
                * while marking this one as a particle node */
            if(firstparent != parent)
            {
                nprnt->f.ChildType = PARTICLE_NODE_TYPE;
                nprnt->s.noccupied = NMAXCHILD;
                tb.Nodes[firstparent].f.ChildType = NODE_NODE_TYPE;
                tb.Nodes[firstparent].s.noccupied = (1<<16);
            }
            return 1;
        }
        for(i=0; i<8; i++) {
            newsuns[i] = newsuns[0] + i;
            struct NODE *nfreep = &tb.Nodes[newsuns[i]];
            /* We create a new leaf node.*/
            init_internal_node(nfreep, nprnt, i);
            /*Set father of new node*/
            nfreep->father = parent;
        }
        /*Initialize the remaining entries to empty*/
        for(i=8; i<NMAXCHILD;i++)
            newsuns[i] = -1;
 
        for(i=0; i < NMAXCHILD; i++) {
            /* Re-attach each particle to the appropriate new leaf.
            * Notice that since we have NMAXCHILD slots on each child and NMAXCHILD particles,
            * we will always have a free slot. */
            int subnode = get_subnode(nprnt, oldsuns[i]);
            int child = newsuns[subnode];
            struct NODE * nchild = &tb.Nodes[child];
            modify_internal_node(child, nchild->s.noccupied, oldsuns[i], tb, HybridNuGrav);
            nchild->s.noccupied++;
        }
        /* Copy the new node array into the node*/
        memcpy(nprnt->s.suns, newsuns, NMAXCHILD * sizeof(int));
        } /* After this brace oldsuns and newsuns are invalid*/
 
        /* Set sibling for the new rank. Since empty at this point, point onwards.*/
        for(i=0; i<7; i++) {
            int child = nprnt->s.suns[i];
            struct NODE * nchild = &tb.Nodes[child];
            nchild->sibling = nprnt->s.suns[i+1];
        }
        /* Final child needs special handling: set to the parent's sibling.*/
        tb.Nodes[nprnt->s.suns[7]].sibling = nprnt->sibling;
        /* Zero the momenta for the parent*/
        memset(&nprnt->mom, 0, sizeof(nprnt->mom));
 
        /* Now try again to add the new particle*/
        int subnode = get_subnode(nprnt, p_toplace);
        int child = nprnt->s.suns[subnode];
        struct NODE * nchild = &tb.Nodes[child];
        if(nchild->s.noccupied < NMAXCHILD) {
            modify_internal_node(child, nchild->s.noccupied, p_toplace, tb, HybridNuGrav);
            nchild->s.noccupied++;
            break;
        }
        /* The attached particles are already within one subnode of the new node.
         * Iterate, creating a new layer beneath.*/
        else {
            /* The current child is going to have new nodes created beneath it,
             * so mark it a Node-containing node. It cannot be accessed until
             * we mark the top-level parent, so no need for atomics.*/
            tb.Nodes[child].f.ChildType = NODE_NODE_TYPE;
            tb.Nodes[child].s.noccupied = (1<<16);
            parent = child;
        }
    } while(1);
 
    /* A new node is created. Mark the (original) parent as an internal node with node children.
     * This goes last so that we don't access the child before it is constructed.*/
    tb.Nodes[firstparent].f.ChildType = NODE_NODE_TYPE;
    tb.Nodes[firstparent].s.noccupied = (1<<16);
    return 0;
}

References NODE::ChildType, NODE::f, NODE::father, get_freenode(), get_subnode(), init_internal_node(), ForceTree::lastnode, message(), modify_internal_node(), NODE::mom, NMAXCHILD, NodeCache::nnext_thread, NodeChild::noccupied, NODE_NODE_TYPE, NODECACHE_SIZE, ForceTree::Nodes, P, PARTICLE_NODE_TYPE, NODE::s, NODE::sibling, and NodeChild::suns.

Referenced by add_particle_to_tree().

Here is the call graph for this function:

Here is the caller graph for this function:

force_create_node_for_topnode()

void force_create_node_for_topnode ( int  no, int  topnode, struct NODE *  Nodes, const DomainDecomp *  ddecomp, int  bits, int  x, int  y, int  z, int *  nextfree, const int  lastnode  )

static

This function recursively creates a set of empty tree nodes which corresponds to the top-level tree for the ddecomp grid. This is done to ensure that this top-level tree is always "complete" so that we can easily associate the pseudo-particles of other CPUs with tree-nodes at a given level in the tree, even when the particle population is so sparse that some of these nodes are actually empty.

Definition at line 818 of file forcetree.c.

{
    int i, j, k;
 
    /*We reached the leaf of the toptree*/
    if(ddecomp->TopNodes[topnode].Daughter < 0)
        return;
 
    for(i = 0; i < 2; i++)
        for(j = 0; j < 2; j++)
            for(k = 0; k < 2; k++)
            {
                int sub = 7 & peano_hilbert_key((x << 1) + i, (y << 1) + j, (z << 1) + k, bits);
 
                int count = i + 2 * j + 4 * k;
 
                Nodes[no].s.Types = 0;
                Nodes[no].s.suns[count] = *nextfree;
                /*We are an internal top level node as we now have a child top level.*/
                Nodes[no].f.InternalTopLevel = 1;
                Nodes[no].f.ChildType = NODE_NODE_TYPE;
                Nodes[no].s.noccupied = (1<<16);
 
                /* We create a new leaf node.*/
                init_internal_node(&Nodes[*nextfree], &Nodes[no], count);
                /*Set father of new node*/
                Nodes[*nextfree].father = no;
                /*All nodes here are top level nodes*/
                Nodes[*nextfree].f.TopLevel = 1;
 
                const struct topnode_data curtopnode = ddecomp->TopNodes[ddecomp->TopNodes[topnode].Daughter + sub];
                if(curtopnode.Daughter == -1) {
                    ddecomp->TopLeaves[curtopnode.Leaf].treenode = *nextfree;
                }
 
                (*nextfree)++;
 
                if(*nextfree >= lastnode)
                    endrun(11, "Not enough force nodes to topnode grid: need %d\n",lastnode);
            }
    /* Set sibling on the child*/
    for(j=0; j<7; j++) {
        int chld = Nodes[no].s.suns[j];
        Nodes[chld].sibling = Nodes[no].s.suns[j+1];
    }
    Nodes[Nodes[no].s.suns[7]].sibling = Nodes[no].sibling;
    for(i = 0; i < 2; i++)
        for(j = 0; j < 2; j++)
            for(k = 0; k < 2; k++)
            {
                int sub = 7 & peano_hilbert_key((x << 1) + i, (y << 1) + j, (z << 1) + k, bits);
                int count = i + 2 * j + 4 * k;
                force_create_node_for_topnode(Nodes[no].s.suns[count], ddecomp->TopNodes[topnode].Daughter + sub, Nodes, ddecomp,
                        bits + 1, 2 * x + i, 2 * y + j, 2 * z + k, nextfree, lastnode);
            }
 
}

References NODE::ChildType, topnode_data::Daughter, endrun(), NODE::f, NODE::father, init_internal_node(), NODE::InternalTopLevel, topnode_data::Leaf, NodeChild::noccupied, NODE_NODE_TYPE, peano_hilbert_key(), NODE::s, NODE::sibling, NodeChild::suns, DomainDecomp::TopLeaves, NODE::TopLevel, DomainDecomp::TopNodes, topleaf_data::treenode, and NodeChild::Types.

Referenced by force_tree_create_nodes().

Here is the call graph for this function:

Here is the caller graph for this function:

force_exchange_pseudodata()

void force_exchange_pseudodata ( ForceTree *  tree, const DomainDecomp *  ddecomp  )

static

This function communicates the values of the multipole moments of the top-level tree-nodes of the ddecomp grid. This data can then be used to update the pseudo-particles on each CPU accordingly.

Definition at line 1109 of file forcetree.c.

{
    int NTask, ThisTask;
    int i, no, ta, recvTask;
    int *recvcounts, *recvoffset;
    struct topleaf_momentsdata
    {
        MyFloat s[3];
        MyFloat mass;
        MyFloat hmax;
    }
    *TopLeafMoments;
 
    MPI_Comm_size(MPI_COMM_WORLD, &NTask);
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
 
    TopLeafMoments = (struct topleaf_momentsdata *) mymalloc("TopLeafMoments", ddecomp->NTopLeaves * sizeof(TopLeafMoments[0]));
    memset(&TopLeafMoments[0], 0, sizeof(TopLeafMoments[0]) * ddecomp->NTopLeaves);
 
    for(i = ddecomp->Tasks[ThisTask].StartLeaf; i < ddecomp->Tasks[ThisTask].EndLeaf; i ++) {
        no = ddecomp->TopLeaves[i].treenode;
        if(ddecomp->TopLeaves[i].Task != ThisTask)
            endrun(131231231, "TopLeaf %d Task table is corrupted: task is %d\n", i, ddecomp->TopLeaves[i].Task);
 
        /* read out the multipole moments from the local base cells */
        TopLeafMoments[i].s[0] = tree->Nodes[no].mom.cofm[0];
        TopLeafMoments[i].s[1] = tree->Nodes[no].mom.cofm[1];
        TopLeafMoments[i].s[2] = tree->Nodes[no].mom.cofm[2];
        TopLeafMoments[i].mass = tree->Nodes[no].mom.mass;
        TopLeafMoments[i].hmax = tree->Nodes[no].mom.hmax;
 
        /*Set the local base nodes dependence on local mass*/
        while(no >= 0)
        {
#ifdef DEBUG
            if(tree->Nodes[no].f.ChildType == PSEUDO_NODE_TYPE)
                endrun(333, "Pseudo node %d parent of a leaf on this processor %d\n", no, ddecomp->TopLeaves[i].treenode);
#endif
            if(tree->Nodes[no].f.DependsOnLocalMass)
                break;
 
            tree->Nodes[no].f.DependsOnLocalMass = 1;
 
            no = tree->Nodes[no].father;
        }
    }
 
    /* share the pseudo-particle data across CPUs */
 
    recvcounts = (int *) mymalloc("recvcounts", sizeof(int) * NTask);
    recvoffset = (int *) mymalloc("recvoffset", sizeof(int) * NTask);
 
    for(recvTask = 0; recvTask < NTask; recvTask++)
    {
        recvoffset[recvTask] = ddecomp->Tasks[recvTask].StartLeaf * sizeof(TopLeafMoments[0]);
        recvcounts[recvTask] = (ddecomp->Tasks[recvTask].EndLeaf - ddecomp->Tasks[recvTask].StartLeaf) * sizeof(TopLeafMoments[0]);
    }
 
    MPI_Allgatherv(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL,
            &TopLeafMoments[0], recvcounts, recvoffset,
            MPI_BYTE, MPI_COMM_WORLD);
 
    myfree(recvoffset);
    myfree(recvcounts);
 
 
    for(ta = 0; ta < NTask; ta++) {
        if(ta == ThisTask) continue; /* bypass ThisTask since it is already up to date */
 
        for(i = ddecomp->Tasks[ta].StartLeaf; i < ddecomp->Tasks[ta].EndLeaf; i ++) {
            no = ddecomp->TopLeaves[i].treenode;
 
            tree->Nodes[no].mom.cofm[0] = TopLeafMoments[i].s[0];
            tree->Nodes[no].mom.cofm[1] = TopLeafMoments[i].s[1];
            tree->Nodes[no].mom.cofm[2] = TopLeafMoments[i].s[2];
            tree->Nodes[no].mom.mass = TopLeafMoments[i].mass;
            tree->Nodes[no].mom.hmax = TopLeafMoments[i].hmax;
         }
    }
    myfree(TopLeafMoments);
}

References NODE::ChildType, NODE::cofm, NODE::DependsOnLocalMass, task_data::EndLeaf, endrun(), NODE::f, NODE::father, NODE::hmax, NODE::mass, NODE::mom, myfree, mymalloc, ForceTree::Nodes, NTask, DomainDecomp::NTopLeaves, PSEUDO_NODE_TYPE, task_data::StartLeaf, topleaf_data::Task, DomainDecomp::Tasks, ThisTask, DomainDecomp::TopLeaves, and topleaf_data::treenode.

Referenced by force_tree_build().

Here is the call graph for this function:

Here is the caller graph for this function:

force_get_father()

int force_get_father	(	int	no,
		const ForceTree *	tree
	)

Definition at line 905 of file forcetree.c.

{
    if(no >= tree->firstnode)
        return tree->Nodes[no].father;
    else
        return tree->Father[no];
}

References NODE::father, ForceTree::Father, ForceTree::firstnode, and ForceTree::Nodes.

Referenced by _prepare(), check_moments(), check_tree(), effhsml(), force_tree_eh_slots_fork(), ngb_treefind_threads(), set_init_hsml(), and setup_smoothinglengths().

Here is the caller graph for this function:

force_get_sibling()

static int force_get_sibling ( const int  sib, const int  j, const int *  suns  )

static

Definition at line 934 of file forcetree.c.

{
    /* check if we have a sibling on the same level */
    int jj;
    int nextsib = sib;
    for(jj = j + 1; jj < 8; jj++) {
        if(suns[jj] >= 0) {
            nextsib = suns[jj];
            break;
        }
    }
    return nextsib;
}

Referenced by force_update_node_recursive().

Here is the caller graph for this function:

force_insert_pseudo_particles()

static void force_insert_pseudo_particles ( const ForceTree *  tree, const DomainDecomp *  ddecomp  )

static

this function inserts pseudo-particles which will represent the mass distribution of the other CPUs. Initially, the mass of the pseudo-particles is set to zero, and their coordinate is set to the center of the ddecomp-cell they correspond to. These quantities will be updated later on.

Definition at line 883 of file forcetree.c.

{
    int i, index;
    const int firstpseudo = tree->lastnode;
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
 
    for(i = 0; i < ddecomp->NTopLeaves; i++)
    {
        index = ddecomp->TopLeaves[i].treenode;
        if(ddecomp->TopLeaves[i].Task != ThisTask) {
            if(tree->Nodes[index].s.noccupied != 0)
                endrun(5, "In node %d, overwriting %d child particles (i = %d etc) with pseudo particle %d\n",
                       index, tree->Nodes[index].s.noccupied, tree->Nodes[index].s.suns[0], i);
            tree->Nodes[index].f.ChildType = PSEUDO_NODE_TYPE;
            /* This node points to the pseudo particle*/
            tree->Nodes[index].s.suns[0] = firstpseudo + i;
        }
    }
}

References NODE::ChildType, endrun(), NODE::f, ForceTree::lastnode, NodeChild::noccupied, ForceTree::Nodes, DomainDecomp::NTopLeaves, PSEUDO_NODE_TYPE, NODE::s, NodeChild::suns, topleaf_data::Task, ThisTask, DomainDecomp::TopLeaves, and topleaf_data::treenode.

Referenced by force_tree_build().

Here is the call graph for this function:

Here is the caller graph for this function:

force_tree_allocated()

int force_tree_allocated

(

const ForceTree *

tree

)

Definition at line 134 of file forcetree.c.

{
    return tree->tree_allocated_flag;
}

References ForceTree::tree_allocated_flag.

Referenced by _prepare(), force_tree_free(), force_tree_rebuild(), force_tree_rebuild_mask(), sfr_reserve_slots(), and treewalk_run().

Here is the caller graph for this function:

force_tree_build()

ForceTree force_tree_build ( int  npart, int  mask, DomainDecomp *  ddecomp, const int  HybridNuGrav, const int  DoMoments, const char *  EmergencyOutputDir  )

static

Constructs the gravitational oct-tree.

The index convention for accessing tree nodes is the following: the indices 0...NumPart-1 reference single particles, the indices firstnode....firstnode +nodes-1 reference tree nodes. ‘Nodes_base’ points to the first tree node, while ‘nodes’ is shifted such that nodes[firstnode] gives the first tree node. Finally, node indices with values 'tb.lastnode' and larger indicate "pseudo particles", i.e. multipole moments of top-level nodes that lie on different CPUs. If such a node needs to be opened, the corresponding particle must be exported to that CPU.

Definition at line 191 of file forcetree.c.

{
    ForceTree tree;
 
    int TooManyNodes = 0;
    int64_t maxnodes = ForceTreeParams.TreeAllocFactor * PartManager->NumPart + ddecomp->NTopNodes;
    int64_t maxmaxnodes;
    MPI_Reduce(&maxnodes, &maxmaxnodes, 1, MPI_INT64, MPI_MAX,0, MPI_COMM_WORLD);
    message(0, "Treebuild: Largest is %g MByte for %ld tree nodes. firstnode %ld. (presently allocated %g MB)\n",
         maxmaxnodes * sizeof(struct NODE) / (1024.0 * 1024.0), maxmaxnodes, PartManager->MaxPart,
         mymalloc_usedbytes() / (1024.0 * 1024.0));
 
    do
    {
        /* Allocate memory. */
        tree = force_treeallocate(maxnodes, PartManager->MaxPart, ddecomp);
 
        tree.BoxSize = PartManager->BoxSize;
        tree.numnodes = force_tree_create_nodes(tree, npart, mask, ddecomp, HybridNuGrav);
        if(tree.numnodes >= tree.lastnode - tree.firstnode)
        {
            message(1, "Not enough tree nodes (%ld) for %d particles. Created %d\n", maxnodes, npart, tree.numnodes);
            force_tree_free(&tree);
            maxnodes = ForceTreeParams.TreeAllocFactor * PartManager->NumPart + ddecomp->NTopNodes;
            message(1, "TreeAllocFactor from %g to %g now %ld tree nodes\n", ForceTreeParams.TreeAllocFactor, ForceTreeParams.TreeAllocFactor*1.15, maxnodes);
            ForceTreeParams.TreeAllocFactor *= 1.15;
            if(ForceTreeParams.TreeAllocFactor > 3.0) {
                TooManyNodes = 1;
                break;
            }
        }
    }
    while(tree.numnodes >= tree.lastnode - tree.firstnode);
 
    tree.mask = mask;
 
    if(MPIU_Any(TooManyNodes, MPI_COMM_WORLD)) {
        /* Assume scale factor = 1 for dump as position is not affected.*/
        if(EmergencyOutputDir) {
            Cosmology CP = {0};
            CP.Omega0 = 0.3;
            CP.OmegaLambda = 0.7;
            CP.HubbleParam = 0.7;
            dump_snapshot("FORCETREE-DUMP", 1, &CP, EmergencyOutputDir);
        }
        endrun(2, "Required too many nodes, snapshot dumped\n");
    }
    if(mask == ALLMASK)
        walltime_measure("/Tree/Build/Nodes");
#ifdef DEBUG
    force_validate_nextlist(&tree);
#endif
    /* insert the pseudo particles that represent the mass distribution of other ddecomps */
    force_insert_pseudo_particles(&tree, ddecomp);
#ifdef DEBUG
    force_validate_nextlist(&tree);
#endif
 
    tree.moments_computed_flag = 0;
 
    if(DoMoments) {
        /* now compute the multipole moments recursively */
        force_update_node_parallel(&tree, ddecomp);
#ifdef DEBUG
        force_validate_nextlist(&tree);
#endif
 
        /* Exchange the pseudo-data*/
        force_exchange_pseudodata(&tree, ddecomp);
 
        force_treeupdate_pseudos(PartManager->MaxPart, &tree);
        tree.moments_computed_flag = 1;
        tree.hmax_computed_flag = 1;
    }
    tree.Nodes_base = (struct NODE *) myrealloc(tree.Nodes_base, (tree.numnodes +1) * sizeof(struct NODE));
 
    /*Update the oct-tree struct so it knows about the memory change*/
    tree.Nodes = tree.Nodes_base - tree.firstnode;
#ifdef DEBUG
        force_validate_nextlist(&tree);
#endif
    return tree;
}

References ALLMASK, ForceTree::BoxSize, part_manager_type::BoxSize, CP, dump_snapshot(), endrun(), ForceTree::firstnode, force_exchange_pseudodata(), force_insert_pseudo_particles(), force_tree_create_nodes(), force_tree_free(), force_treeallocate(), force_treeupdate_pseudos(), force_update_node_parallel(), ForceTreeParams, ForceTree::hmax_computed_flag, Cosmology::HubbleParam, ForceTree::lastnode, ForceTree::mask, part_manager_type::MaxPart, message(), ForceTree::moments_computed_flag, MPI_INT64, MPIU_Any(), mymalloc_usedbytes, myrealloc, ForceTree::Nodes, ForceTree::Nodes_base, DomainDecomp::NTopNodes, ForceTree::numnodes, part_manager_type::NumPart, Cosmology::Omega0, Cosmology::OmegaLambda, PartManager, forcetree_params::TreeAllocFactor, and walltime_measure.

Referenced by force_tree_rebuild(), and force_tree_rebuild_mask().

Here is the call graph for this function:

Here is the caller graph for this function:

force_tree_create_nodes()

int force_tree_create_nodes	(	const ForceTree	tb,
		const int	npart,
		int	mask,
		DomainDecomp *	ddecomp,
		const int	HybridNuGrav
	)

Does initial creation of the nodes for the gravitational oct-tree. mask is a bitfield: Only types whose bit is set are added.

Definition at line 661 of file forcetree.c.

{
    int nnext = tb.firstnode;       /* index of first free node */
 
    /* create an empty root node  */
    {
        int i;
        struct NODE *nfreep = &tb.Nodes[nnext]; /* select first node */
 
        nfreep->len = PartManager->BoxSize*1.001;
        for(i = 0; i < 3; i++)
            nfreep->center[i] = PartManager->BoxSize/2.;
        for(i = 0; i < NMAXCHILD; i++)
            nfreep->s.suns[i] = -1;
        nfreep->s.Types = 0;
        nfreep->s.noccupied = 0;
        nfreep->father = -1;
        nfreep->sibling = -1;
        nfreep->f.TopLevel = 1;
        nfreep->f.InternalTopLevel = 0;
        nfreep->f.DependsOnLocalMass = 0;
        nfreep->f.ChildType = PARTICLE_NODE_TYPE;
        nfreep->f.unused = 0;
        memset(&(nfreep->mom.cofm),0,3*sizeof(MyFloat));
        nfreep->mom.mass = 0;
        nfreep->mom.hmax = 0;
        nnext++;
        /* create a set of empty nodes corresponding to the top-level ddecomp
         * grid. We need to generate these nodes first to make sure that we have a
         * complete top-level tree which allows the easy insertion of the
         * pseudo-particles in the right place */
        force_create_node_for_topnode(tb.firstnode, 0, tb.Nodes, ddecomp, 1, 0, 0, 0, &nnext, tb.lastnode);
    }
    /* Set up thread-local copies of the topnodes to anchor the subtrees. */
    int ThisTask, j, t;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    const int StartLeaf = ddecomp->Tasks[ThisTask].StartLeaf;
    const int EndLeaf = ddecomp->Tasks[ThisTask].EndLeaf;
    const int nthr = omp_get_max_threads();
    int * topnodes = ta_malloc("topnodes", int, (EndLeaf - StartLeaf) * nthr);
    /* Topnodes for each thread. For tid 0, just use the real tree. Saves copying the tree back.*/
    for(j = 0; j < EndLeaf - StartLeaf; j++)
        topnodes[j] = ddecomp->TopLeaves[j + StartLeaf].treenode;
    /* Other threads need a copy*/
    for(t = 1; t < nthr; t++) {
        for(j = 0; j < EndLeaf - StartLeaf; j++) {
            /* Make a local copy*/
            topnodes[j + t * (EndLeaf - StartLeaf)] = nnext;
            memmove(&tb.Nodes[nnext], &tb.Nodes[topnodes[j]], sizeof(struct NODE));
            nnext++;
        }
    }
 
/*     double tstart = second(); */
    const int first_free = nnext;
    /* Increment nnext for the threads we are about to initialise.*/
    nnext += NODECACHE_SIZE * nthr;
    /* now we insert all particles */
    #pragma omp parallel
    {
        int i;
        /* Local topnodes*/
        int tid = omp_get_thread_num();
        const int * const local_topnodes = topnodes + tid * (EndLeaf - StartLeaf);
 
        /* This implements a small thread-local free Node cache.
         * The cache ensures that Nodes from the same (or close) particles
         * are created close to each other on the Node list and thus
         * helps cache locality. I tried each thread getting a separate
         * part of the tree, and it wasted too much memory. */
        struct NodeCache nc;
        nc.nnext_thread = first_free + tid * NODECACHE_SIZE;
        nc.nrem_thread = NODECACHE_SIZE;
 
        /* Stores the last-seen node on this thread.
         * Since most particles are close to each other, this should save a number of tree walks.*/
        int this_acc = local_topnodes[0];
        // message(1, "Topnodes %d real %d\n", local_topnodes[0], topnodes[0]);
 
        /* The default schedule is static with a chunk 1/4 the total.
         * However, particles are sorted by type and then by peano order.
         * Since we need to merge trees, it is advantageous to have all particles
         * spatially close be processed by the same thread. This means that the threads should
         * process particles at a constant offset from the start of the type.
         * We do this with a static schedule. */
        int chnksz = PartManager->NumPart/nthr;
        if(SlotsManager->info[0].enabled && SlotsManager->info[0].size > 0)
            chnksz = SlotsManager->info[0].size/nthr;
        if(chnksz < 1000)
            chnksz = 1000;
        #pragma omp for schedule(static, chnksz)
        for(i = 0; i < npart; i++)
        {
            /*Can't break from openmp for*/
            if(nc.nnext_thread >= tb.lastnode)
                continue;
 
            /* Do not add types that do not have their mask bit set.*/
            if(!((1<<P[i].Type) & mask)) {
                continue;
            }
            /* Do not add garbage/swallowed particles to the tree*/
            if(P[i].IsGarbage || (P[i].Swallowed && P[i].Type==5))
                continue;
 
            if(P[i].Mass == 0)
                endrun(12, "Zero mass particle %d type %d id %ld pos %g %g %g\n", i, P[i].Type, P[i].ID, P[i].Pos[0], P[i].Pos[1], P[i].Pos[2]);
            /*First find the Node for the TopLeaf */
            int cur;
            if(inside_node(&tb.Nodes[this_acc], i)) {
                cur = this_acc;
            } else {
                /* Get the topnode to which a particle belongs. Each local tree
                 * has a local set of treenodes copying the global topnodes, except tid 0
                 * which has the real topnodes.*/
                const int topleaf = domain_get_topleaf(P[i].Key, ddecomp);
                //int treenode = ddecomp->TopLeaves[topleaf].treenode;
                cur = local_topnodes[topleaf - StartLeaf];
            }
 
            this_acc = add_particle_to_tree(i, cur, tb, HybridNuGrav, &nc, &nnext);
        }
        /* The implicit omp-barrier is important here!*/
/*         double tend = second(); */
/*         message(0, "Initial insertion: %.3g ms. First node %d\n", (tend - tstart)*1000, local_topnodes[0]); */
 
        /* Merge each topnode separately, using a for loop.
         * This wastes threads if NTHREAD > NTOPNODES, but it
         * means only one merge is done per subtree and
         * it requires no locking.*/
        #pragma omp for schedule(static, 1)
        for(i = 0; i < EndLeaf - StartLeaf; i++) {
            int t;
            /* These are the addresses of the real topnodes*/
            const int target = topnodes[i];
            if(nc.nnext_thread >= tb.lastnode)
                continue;
            for(t = 1; t < nthr; t++) {
                const int righttop = topnodes[i + t * (EndLeaf - StartLeaf)];
//                  message(1, "tid = %d i = %d t = %d Merging %d to %d addresses are %lx - %lx end is %lx\n", omp_get_thread_num(), i, t, righttop, target, &tb.Nodes[righttop], &tb.Nodes[target], &tb.Nodes[nnext]);
                if(merge_partial_force_trees(target, righttop, &nc, &nnext, tb, HybridNuGrav))
                    break;
            }
        }
    }
 
    ta_free(topnodes);
    return nnext - tb.firstnode;
}

References add_particle_to_tree(), part_manager_type::BoxSize, NODE::center, NODE::ChildType, NODE::cofm, NODE::DependsOnLocalMass, domain_get_topleaf(), slot_info::enabled, task_data::EndLeaf, endrun(), NODE::f, NODE::father, ForceTree::firstnode, force_create_node_for_topnode(), NODE::hmax, slots_manager_type::info, inside_node(), NODE::InternalTopLevel, ForceTree::lastnode, NODE::len, NODE::mass, merge_partial_force_trees(), NODE::mom, NMAXCHILD, NodeCache::nnext_thread, NodeChild::noccupied, NODECACHE_SIZE, ForceTree::Nodes, NodeCache::nrem_thread, part_manager_type::NumPart, P, PARTICLE_NODE_TYPE, PartManager, NODE::s, NODE::sibling, slot_info::size, SlotsManager, task_data::StartLeaf, NodeChild::suns, ta_free, ta_malloc, DomainDecomp::Tasks, ThisTask, DomainDecomp::TopLeaves, NODE::TopLevel, topleaf_data::treenode, NodeChild::Types, and NODE::unused.

Referenced by do_tree_test(), and force_tree_build().

Here is the call graph for this function:

Here is the caller graph for this function:

force_tree_eh_slots_fork()

static int force_tree_eh_slots_fork ( EIBase *  event, void *  userdata  )

static

Definition at line 110 of file forcetree.c.

{
    /* after a fork, we will attach the new particle to the force tree. */
    EISlotsFork * ev = (EISlotsFork*) event;
    int parent = ev->parent;
    int child = ev->child;
    ForceTree * tree = (ForceTree * ) userdata;
    int no = force_get_father(parent, tree);
    struct NODE * nop = &tree->Nodes[no];
    /* FIXME: We lose particles if the node is full.
     * At the moment this does not matter, because
     * the only new particles are stars, which do not
     * participate in the SPH tree walk.*/
    if(nop->s.noccupied < NMAXCHILD) {
       nop->s.suns[nop->s.noccupied] = child;
       nop->s.Types += P[child].Type << (3*nop->s.noccupied);
       nop->s.noccupied++;
    }
    if(child < tree->nfather && child >= 0)
        tree->Father[child] = no;
    return 0;
}

References EISlotsFork::child, ForceTree::Father, force_get_father(), NMAXCHILD, NodeChild::noccupied, ForceTree::Nodes, P, EISlotsFork::parent, NODE::s, NodeChild::suns, and NodeChild::Types.

Referenced by force_tree_free(), and force_tree_rebuild().

Here is the call graph for this function:

Here is the caller graph for this function:

force_tree_free()

void force_tree_free

(

ForceTree *

tree

)

This function frees the memory allocated for the tree, i.e. it frees the space allocated by the function force_treeallocate().

Definition at line 1409 of file forcetree.c.

{
    event_unlisten(&EventSlotsFork, force_tree_eh_slots_fork, tree);
 
    if(!force_tree_allocated(tree))
        return;
    myfree(tree->Nodes_base);
    myfree(tree->Father);
    tree->tree_allocated_flag = 0;
}

References event_unlisten(), EventSlotsFork, ForceTree::Father, force_tree_allocated(), force_tree_eh_slots_fork(), myfree, ForceTree::Nodes_base, and ForceTree::tree_allocated_flag.

Referenced by _prepare(), do_density_test(), do_force_test(), fof_fof(), force_tree_build(), force_tree_rebuild(), force_tree_rebuild_mask(), metal_return(), run(), run_gravity_test(), runfof(), runpower(), setup_smoothinglengths(), test_rebuild_close(), test_rebuild_flat(), test_rebuild_random(), and turn_on_quasars().

Here is the call graph for this function:

Here is the caller graph for this function:

force_tree_rebuild()

void force_tree_rebuild	(	ForceTree *	tree,
		DomainDecomp *	ddecomp,
		const int	HybridNuGrav,
		const int	DoMoments,
		const char *	EmergencyOutputDir
	)

Definition at line 140 of file forcetree.c.

{
    MPIU_Barrier(MPI_COMM_WORLD);
    message(0, "Tree construction.  (presently allocated=%g MB)\n", mymalloc_usedbytes() / (1024.0 * 1024.0));
 
    if(force_tree_allocated(tree)) {
        force_tree_free(tree);
    }
    walltime_measure("/Misc");
 
    *tree = force_tree_build(PartManager->NumPart, ALLMASK, ddecomp, HybridNuGrav, DoMoments, EmergencyOutputDir);
 
    event_listen(&EventSlotsFork, force_tree_eh_slots_fork, tree);
    walltime_measure("/Tree/Build/Moments");
 
    message(0, "Tree constructed (moments: %d). First node %d, number of nodes %d, first pseudo %d. NTopLeaves %d\n",
            tree->moments_computed_flag, tree->firstnode, tree->numnodes, tree->lastnode, tree->NTopLeaves);
    MPIU_Barrier(MPI_COMM_WORLD);
}

References ALLMASK, event_listen(), EventSlotsFork, ForceTree::firstnode, force_tree_allocated(), force_tree_build(), force_tree_eh_slots_fork(), force_tree_free(), ForceTree::lastnode, message(), ForceTree::moments_computed_flag, MPIU_Barrier, mymalloc_usedbytes, ForceTree::NTopLeaves, ForceTree::numnodes, part_manager_type::NumPart, PartManager, and walltime_measure.

Referenced by do_density_test(), do_force_test(), run(), run_gravity_test(), runfof(), runpower(), and setup_smoothinglengths().

Here is the call graph for this function:

Here is the caller graph for this function:

force_tree_rebuild_mask()

void force_tree_rebuild_mask	(	ForceTree *	tree,
		DomainDecomp *	ddecomp,
		int	mask,
		const int	HybridNuGrav,
		const char *	EmergencyOutputDir
	)

Definition at line 161 of file forcetree.c.

{
    MPIU_Barrier(MPI_COMM_WORLD);
    message(0, "Tree construction for types: %d.\n", mask);
 
    if(force_tree_allocated(tree)) {
        force_tree_free(tree);
    }
    walltime_measure("/Misc");
 
    /* No moments*/
    *tree = force_tree_build(PartManager->NumPart, mask, ddecomp, HybridNuGrav, 0, EmergencyOutputDir);
    walltime_measure("/SPH/Build");
 
    message(0, "Tree constructed (type mask: %d). First node %d, number of nodes %d, first pseudo %d. NTopLeaves %d\n",
            mask, tree->firstnode, tree->numnodes, tree->lastnode, tree->NTopLeaves);
    MPIU_Barrier(MPI_COMM_WORLD);
}

References ForceTree::firstnode, force_tree_allocated(), force_tree_build(), force_tree_free(), ForceTree::lastnode, message(), MPIU_Barrier, ForceTree::NTopLeaves, ForceTree::numnodes, part_manager_type::NumPart, PartManager, and walltime_measure.

Referenced by fof_fof(), metal_return(), and turn_on_quasars().

Here is the call graph for this function:

Here is the caller graph for this function:

force_treeallocate()

ForceTree force_treeallocate	(	int64_t	maxnodes,
		int64_t	maxpart,
		DomainDecomp *	ddecomp
	)

This function allocates the memory used for storage of the tree and of auxiliary arrays needed for tree-walk and link-lists. Usually, maxnodes approximately equal to 0.7*maxpart is sufficient to store the tree for up to maxpart particles.

Definition at line 1381 of file forcetree.c.

{
    ForceTree tb;
 
    tb.Father = (int *) mymalloc("Father", maxpart * sizeof(int));
    tb.nfather = maxpart;
#ifdef DEBUG
    memset(tb.Father, -1, maxpart * sizeof(int));
#endif
    tb.Nodes_base = (struct NODE *) mymalloc("Nodes_base", (maxnodes + 1) * sizeof(struct NODE));
#ifdef DEBUG
    memset(tb.Nodes_base, -1, (maxnodes + 1) * sizeof(struct NODE));
#endif
    tb.firstnode = maxpart;
    tb.lastnode = maxpart + maxnodes;
    if(maxpart + maxnodes >= 1L<<30)
        endrun(5, "Size of tree overflowed for maxpart = %ld, maxnodes = %ld!\n", maxpart, maxnodes);
    tb.numnodes = 0;
    tb.Nodes = tb.Nodes_base - maxpart;
    tb.tree_allocated_flag = 1;
    tb.NTopLeaves = ddecomp->NTopLeaves;
    tb.TopLeaves = ddecomp->TopLeaves;
    return tb;
}

References endrun(), ForceTree::Father, ForceTree::firstnode, ForceTree::lastnode, mymalloc, ForceTree::nfather, ForceTree::Nodes, ForceTree::Nodes_base, DomainDecomp::NTopLeaves, ForceTree::NTopLeaves, ForceTree::numnodes, DomainDecomp::TopLeaves, ForceTree::TopLeaves, and ForceTree::tree_allocated_flag.

Referenced by force_tree_build(), test_rebuild_close(), test_rebuild_flat(), and test_rebuild_random().

Here is the call graph for this function:

Here is the caller graph for this function:

force_treeupdate_pseudos()

void force_treeupdate_pseudos ( int  no, const ForceTree *  tree  )

static

This function updates the top-level tree after the multipole moments of the pseudo-particles have been updated.

Definition at line 1195 of file forcetree.c.

{
    int j, p;
    MyFloat hmax;
    MyFloat s[3], mass;
 
    mass = 0;
    s[0] = 0;
    s[1] = 0;
    s[2] = 0;
    hmax = 0;
 
    /* This happens if we have a trivial domain with only one entry*/
    if(!tree->Nodes[no].f.InternalTopLevel)
        return;
 
    p = tree->Nodes[no].s.suns[0];
 
    /* since we are dealing with top-level nodes, we know that there are 8 consecutive daughter nodes */
    for(j = 0; j < 8; j++)
    {
        /*This may not happen as we are an internal top level node*/
        if(p < tree->firstnode || p >= tree->lastnode)
            endrun(6767, "Updating pseudos: %d -> %d which is not an internal node between %d and %d.",no, p, tree->firstnode, tree->lastnode);
#ifdef DEBUG
        /* Check we don't move to another part of the tree*/
        if(tree->Nodes[p].father != no)
            endrun(6767, "Tried to update toplevel node %d with parent %d != expected %d\n", p, tree->Nodes[p].father, no);
#endif
 
        if(tree->Nodes[p].f.InternalTopLevel)
            force_treeupdate_pseudos(p, tree);
 
        mass += (tree->Nodes[p].mom.mass);
        s[0] += (tree->Nodes[p].mom.mass * tree->Nodes[p].mom.cofm[0]);
        s[1] += (tree->Nodes[p].mom.mass * tree->Nodes[p].mom.cofm[1]);
        s[2] += (tree->Nodes[p].mom.mass * tree->Nodes[p].mom.cofm[2]);
 
        if(tree->Nodes[p].mom.hmax > hmax)
            hmax = tree->Nodes[p].mom.hmax;
 
        p = tree->Nodes[p].sibling;
    }
 
    if(mass)
    {
        s[0] /= mass;
        s[1] /= mass;
        s[2] /= mass;
    }
    else
    {
        s[0] = tree->Nodes[no].center[0];
        s[1] = tree->Nodes[no].center[1];
        s[2] = tree->Nodes[no].center[2];
    }
 
    tree->Nodes[no].mom.cofm[0] = s[0];
    tree->Nodes[no].mom.cofm[1] = s[1];
    tree->Nodes[no].mom.cofm[2] = s[2];
    tree->Nodes[no].mom.mass = mass;
 
    tree->Nodes[no].mom.hmax = hmax;
}

References NODE::center, NODE::cofm, endrun(), NODE::f, NODE::father, ForceTree::firstnode, NODE::hmax, NODE::InternalTopLevel, ForceTree::lastnode, NODE::mass, NODE::mom, ForceTree::Nodes, NODE::s, NODE::sibling, and NodeChild::suns.

Referenced by force_tree_build().

Here is the call graph for this function:

Here is the caller graph for this function:

force_update_hmax()

void force_update_hmax	(	int *	activeset,
		int	size,
		ForceTree *	tree,
		DomainDecomp *	ddecomp
	)

This function updates the hmax-values in tree nodes that hold SPH particles. Since the Hsml-values are potentially changed for active particles in the SPH-density computation, force_update_hmax() should be carried out just before the hydrodynamical SPH forces are computed, i.e. after density().

The purpose of the hmax node is for a symmetric treewalk (currently only the hydro). Particles where P[i].Pos + Hsml pokes beyond the exterior of the tree node may mean that a tree node should be included when it would normally be culled. Therefore we don't really want hmax, we want the maximum amount P[i].Pos + Hsml pokes beyond the tree node.

Definition at line 1271 of file forcetree.c.

{
    int NTask, ThisTask, recvTask;
    int i, ta;
    int *recvcounts, *recvoffset;
 
    walltime_measure("/Misc");
 
    /* If hmax has not yet been computed, do all particles*/
    if(!tree->hmax_computed_flag) {
        activeset = NULL;
        size = PartManager->NumPart;
    }
    MPI_Comm_size(MPI_COMM_WORLD, &NTask);
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
 
    #pragma omp parallel for
    for(i = 0; i < size; i++)
    {
        const int p_i = activeset ? activeset[i] : i;
 
        if(P[p_i].Type != 0 || P[p_i].IsGarbage)
            continue;
 
        int no = tree->Father[p_i];
 
        while(no >= 0)
        {
            /* How much does this particle peek beyond this node?
             * Note len does not change so we can read it without a lock or atomic. */
            MyFloat readhmax, newhmax = 0;
            int j, done = 0;
            for(j = 0; j < 3; j++) {
                /* Compute each direction independently and take the maximum.
                 * This is the largest possible distance away from node center within a cube bounding hsml.
                 * Note that because Pos - Center < len, the maximum value this can have is Hsml.*/
                newhmax = DMAX(newhmax, fabs(P[p_i].Pos[j] - tree->Nodes[no].center[j]) + P[p_i].Hsml - tree->Nodes[no].len);
            }
            /* Most particles will lie fully inside a node. No need then for the atomic! */
            if(newhmax <= 0)
                break;
 
            #pragma omp atomic read
            readhmax = tree->Nodes[no].mom.hmax;
            do {
                if(newhmax <= readhmax) {
                    done = 1;
                    break;
                }
                /* Swap in the new hmax only if the old one hasn't changed. */
            } while(!__atomic_compare_exchange(&(tree->Nodes[no].mom.hmax), &readhmax, &newhmax, 0, __ATOMIC_RELAXED, __ATOMIC_RELAXED));
 
            if(done)
                break;
            no = tree->Nodes[no].father;
        }
    }
 
    double * TopLeafhmax = (double *) mymalloc("TopLeafMoments", ddecomp->NTopLeaves * sizeof(double));
    memset(&TopLeafhmax[0], 0, sizeof(double) * ddecomp->NTopLeaves);
 
    for(i = ddecomp->Tasks[ThisTask].StartLeaf; i < ddecomp->Tasks[ThisTask].EndLeaf; i ++) {
        int no = ddecomp->TopLeaves[i].treenode;
        TopLeafhmax[i] = tree->Nodes[no].mom.hmax;
    }
 
    /* share the hmax-data of the dirty nodes accross CPUs */
    recvcounts = (int *) mymalloc("recvcounts", sizeof(int) * NTask);
    recvoffset = (int *) mymalloc("recvoffset", sizeof(int) * NTask);
 
    for(recvTask = 0; recvTask < NTask; recvTask++)
    {
        recvoffset[recvTask] = ddecomp->Tasks[recvTask].StartLeaf;
        recvcounts[recvTask] = ddecomp->Tasks[recvTask].EndLeaf - ddecomp->Tasks[recvTask].StartLeaf;
    }
 
    MPI_Allgatherv(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL,
            &TopLeafhmax[0], recvcounts, recvoffset,
            MPI_DOUBLE, MPI_COMM_WORLD);
 
    myfree(recvoffset);
    myfree(recvcounts);
 
    for(ta = 0; ta < NTask; ta++) {
        if(ta == ThisTask)
            continue; /* bypass ThisTask since it is already up to date */
        for(i = ddecomp->Tasks[ta].StartLeaf; i < ddecomp->Tasks[ta].EndLeaf; i ++) {
            int no = ddecomp->TopLeaves[i].treenode;
            tree->Nodes[no].mom.hmax = TopLeafhmax[i];
 
            while(no >= 0)
            {
                if(TopLeafhmax[i] <= tree->Nodes[no].mom.hmax)
                    break;
                tree->Nodes[no].mom.hmax = TopLeafhmax[i];
                no = tree->Nodes[no].father;
            }
         }
    }
    myfree(TopLeafhmax);
 
    tree->hmax_computed_flag = 1;
    walltime_measure("/Tree/HmaxUpdate");
}

References NODE::center, DMAX, task_data::EndLeaf, NODE::father, ForceTree::Father, NODE::hmax, ForceTree::hmax_computed_flag, NODE::len, NODE::mom, myfree, mymalloc, ForceTree::Nodes, NTask, DomainDecomp::NTopLeaves, part_manager_type::NumPart, P, PartManager, task_data::StartLeaf, DomainDecomp::Tasks, ThisTask, DomainDecomp::TopLeaves, topleaf_data::treenode, and walltime_measure.

Referenced by run().

Here is the caller graph for this function:

force_update_node_parallel()

void force_update_node_parallel	(	const ForceTree *	tree,
		const DomainDecomp *	ddecomp
	)

This routine determines the multipole moments for a given internal node and all its subnodes in parallel, assigning the recursive algorithm to different threads using openmp's task api. The result is stored in tb.Nodes in the sequence of this tree-walk.

• A new task is spawned from each down from each local topleaf. Local topleaves are used so that we do not waste time trying moment calculation with pseudoparticles.
• Each internal node found at that level is added to a list, together with its sibling.
• Each node in this list then has the recursive moment calculation called on it. Note: If the tree is very unbalanced and one branch much deeper than the others, this will not be efficient.
• Once each tree's recursive moment is generated in parallel, the tail value from each recursion is stored, and the node marked as done.
• A final recursive moment calculation is run in serial for the top 3 levels of the tree. When it encounters one of the pre-computed nodes, it searches the list of pre-computed tail values to set the next node as if it had recursed and continues.

Definition at line 1081 of file forcetree.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
 
#pragma omp parallel
#pragma omp single nowait
    {
        int i;
        for(i = ddecomp->Tasks[ThisTask].StartLeaf; i < ddecomp->Tasks[ThisTask].EndLeaf; i ++) {
            const int no = ddecomp->TopLeaves[i].treenode;
            /* Nodes containing other nodes: the overwhelmingly likely case.*/
            if(tree->Nodes[no].f.ChildType == NODE_NODE_TYPE) {
                #pragma omp task default(none) shared(tree) firstprivate(no)
                force_update_node_recursive(no, tree->Nodes[no].sibling, 1, tree);
            }
            else if(tree->Nodes[no].f.ChildType == PARTICLE_NODE_TYPE)
                force_update_particle_node(no, tree);
            else if(tree->Nodes[no].f.ChildType == PSEUDO_NODE_TYPE)
                endrun(5, "Error, found pseudo node %d but domain entry %d says on task %d\n", no, i, ThisTask);
        }
    }
}

References NODE::ChildType, endrun(), NODE::f, force_update_node_recursive(), force_update_particle_node(), NODE_NODE_TYPE, ForceTree::Nodes, PARTICLE_NODE_TYPE, PSEUDO_NODE_TYPE, NODE::sibling, task_data::StartLeaf, DomainDecomp::Tasks, ThisTask, DomainDecomp::TopLeaves, and topleaf_data::treenode.

Referenced by do_tree_test(), and force_tree_build().

Here is the call graph for this function:

Here is the caller graph for this function:

force_update_node_recursive()

static int force_update_node_recursive ( int  no, int  sib, int  level, const ForceTree *  tree  )

static

this routine determines the multipole moments for a given internal node and all its subnodes using a recursive computation. The result is stored in tb.Nodes in the sequence of this tree-walk.

The function also computes the NextNode and sibling linked lists. The return value is the current tail of the NextNode linked list.

This function is called recursively using openmp tasks. We spawn a new task for a fixed number of levels of the tree.

Definition at line 982 of file forcetree.c.

{
#ifdef DEBUG
    if(tree->Nodes[no].f.ChildType != NODE_NODE_TYPE)
        endrun(3, "force_update_node_recursive called on node %d of type %d != %d!\n", no, tree->Nodes[no].f.ChildType, NODE_NODE_TYPE);
#endif
    int j;
    int * suns = tree->Nodes[no].s.suns;
 
    int childcnt = 0;
    /* Remove any empty children, moving the suns array around
     * so non-empty entries are contiguous at the beginning of the array.
     * This sharply reduces the size of the tree.
     * Also count the node children for thread balancing.*/
    int jj = 0;
    for(j=0; j < 8; j++, jj++) {
        /* Never remove empty top-level nodes so we don't
         * mess up the pseudo-data exchange.
         * This may happen for a pseudo particle host or, in very rare cases,
         * when one of the local domains is empty. */
        while(jj < 8 && !tree->Nodes[suns[jj]].f.TopLevel &&
            tree->Nodes[suns[jj]].f.ChildType == PARTICLE_NODE_TYPE &&
            tree->Nodes[suns[jj]].s.noccupied == 0) {
                    jj++;
        }
        if(jj < 8)
            suns[j] = suns[jj];
        else
            suns[j] = -1;
        if(suns[j] >= 0 && tree->Nodes[suns[j]].f.ChildType == NODE_NODE_TYPE)
            childcnt++;
    }
 
    /*First do the children*/
    for(j = 0; j < 8; j++)
    {
        int p = suns[j];
        /*Empty slot*/
        if(p < 0)
            continue;
        const int nextsib = force_get_sibling(sib, j, suns);
        /* This is set in create_nodes but needed because we may remove empty nodes above.*/
        tree->Nodes[p].sibling = nextsib;
        /* Nodes containing particles or pseudo-particles*/
        if(tree->Nodes[p].f.ChildType == PARTICLE_NODE_TYPE)
            force_update_particle_node(p, tree);
        if(tree->Nodes[p].f.ChildType == NODE_NODE_TYPE) {
            /* Don't spawn a new task if we are deep enough that we already spawned a lot.*/
            if(childcnt > 1 && level < 512) {
                #pragma omp task default(none) shared(level, childcnt, tree) firstprivate(nextsib, p)
                force_update_node_recursive(p, nextsib, level*childcnt, tree);
            }
            else
                force_update_node_recursive(p, nextsib, level, tree);
        }
    }
 
    /*Make sure all child nodes are done*/
    #pragma omp taskwait
 
    /*Now we do the moments*/
    for(j = 0; j < 8; j++)
    {
        const int p = suns[j];
        if(p < 0)
            continue;
        tree->Nodes[no].mom.mass += (tree->Nodes[p].mom.mass);
        tree->Nodes[no].mom.cofm[0] += (tree->Nodes[p].mom.mass * tree->Nodes[p].mom.cofm[0]);
        tree->Nodes[no].mom.cofm[1] += (tree->Nodes[p].mom.mass * tree->Nodes[p].mom.cofm[1]);
        tree->Nodes[no].mom.cofm[2] += (tree->Nodes[p].mom.mass * tree->Nodes[p].mom.cofm[2]);
        if(tree->Nodes[p].mom.hmax > tree->Nodes[no].mom.hmax)
            tree->Nodes[no].mom.hmax = tree->Nodes[p].mom.hmax;
    }
 
    /*Set the center of mass moments*/
    const double mass = tree->Nodes[no].mom.mass;
    /* In principle all the children could be pseudo-particles*/
    if(mass > 0) {
        tree->Nodes[no].mom.cofm[0] /= mass;
        tree->Nodes[no].mom.cofm[1] /= mass;
        tree->Nodes[no].mom.cofm[2] /= mass;
    }
 
    return -1;
}

References NODE::ChildType, NODE::cofm, endrun(), NODE::f, force_get_sibling(), force_update_particle_node(), NODE::hmax, NODE::mass, NODE::mom, NodeChild::noccupied, NODE_NODE_TYPE, ForceTree::Nodes, PARTICLE_NODE_TYPE, NODE::s, NODE::sibling, and NodeChild::suns.

Referenced by force_update_node_parallel().

Here is the call graph for this function:

Here is the caller graph for this function:

force_update_particle_node()

static void force_update_particle_node ( int  no, const ForceTree *  tree  )

static

Definition at line 950 of file forcetree.c.

{
#ifdef DEBUG
    if(tree->Nodes[no].f.ChildType != PARTICLE_NODE_TYPE)
        endrun(3, "force_update_particle_node called on node %d of wrong type!\n", no);
#endif
    int j;
    /*Set the center of mass moments*/
    const double mass = tree->Nodes[no].mom.mass;
    /* Be careful about empty nodes*/
    if(mass > 0) {
        for(j = 0; j < 3; j++)
            tree->Nodes[no].mom.cofm[j] /= mass;
    }
    else {
        for(j = 0; j < 3; j++)
            tree->Nodes[no].mom.cofm[j] = tree->Nodes[no].center[j];
    }
}

References NODE::center, NODE::ChildType, NODE::cofm, endrun(), NODE::f, NODE::mass, NODE::mom, ForceTree::Nodes, and PARTICLE_NODE_TYPE.

Referenced by force_update_node_parallel(), and force_update_node_recursive().

Here is the call graph for this function:

Here is the caller graph for this function:

get_freenode()

int get_freenode	(	int *	nnext,
		struct NodeCache *	nc
	)

Definition at line 345 of file forcetree.c.

{
    /*Get memory for an extra node from our cache.*/
    if(nc->nrem_thread < 8) {
        nc->nnext_thread = atomic_fetch_and_add(nnext, NODECACHE_SIZE);
        nc->nrem_thread = NODECACHE_SIZE;
    }
    const int ninsert = nc->nnext_thread;
    nc->nnext_thread += 8;
    nc->nrem_thread -= 8;
    return ninsert;
}

References atomic_fetch_and_add(), NodeCache::nnext_thread, NODECACHE_SIZE, and NodeCache::nrem_thread.

Referenced by create_new_node_layer().

Here is the call graph for this function:

Here is the caller graph for this function:

get_subnode()

int get_subnode	(	const struct NODE *	node,
		const int	p_i
	)

Definition at line 281 of file forcetree.c.

{
    /*Loop is unrolled to help out the compiler,which normally only manages it at -O3*/
     return (P[p_i].Pos[0] > node->center[0]) +
            ((P[p_i].Pos[1] > node->center[1]) << 1) +
            ((P[p_i].Pos[2] > node->center[2]) << 2);
}

References NODE::center, and P.

Referenced by add_particle_to_tree(), and create_new_node_layer().

Here is the caller graph for this function:

init_forcetree_params()

void init_forcetree_params

(

const int

FastParticleType

)

Definition at line 43 of file forcetree.c.

{
    /* This was increased due to the extra nodes created by subtrees*/
    ForceTreeParams.TreeAllocFactor = 0.9;
    ForceTreeParams.FastParticleType = FastParticleType;
}

References forcetree_params::FastParticleType, ForceTreeParams, and forcetree_params::TreeAllocFactor.

Referenced by begrun(), setup_density(), setup_tree(), and test_fof().

Here is the caller graph for this function:

init_internal_node()

static void init_internal_node ( struct NODE *  nfreep, struct NODE *  parent, int  subnode  )

static

Definition at line 304 of file forcetree.c.

{
    int j;
    const MyFloat lenhalf = 0.25 * parent->len;
    nfreep->len = 0.5 * parent->len;
    nfreep->sibling = -10;
    nfreep->father = -10;
    nfreep->f.TopLevel = 0;
    nfreep->f.InternalTopLevel = 0;
    nfreep->f.DependsOnLocalMass = 0;
    nfreep->f.ChildType = PARTICLE_NODE_TYPE;
    nfreep->f.unused = 0;
 
    for(j = 0; j < 3; j++) {
        /* Detect which quadrant we are in by testing the bits of subnode:
         * if (subnode & [1,2,4]) is true we add lenhalf, otherwise subtract lenhalf*/
        const int sign = (subnode & (1 << j)) ? 1 : -1;
        nfreep->center[j] = parent->center[j] + sign*lenhalf;
    }
    for(j = 0; j < NMAXCHILD; j++)
        nfreep->s.suns[j] = -1;
    nfreep->s.noccupied = 0;
    nfreep->s.Types = 0;
    memset(&(nfreep->mom.cofm),0,3*sizeof(MyFloat));
    nfreep->mom.mass = 0;
    nfreep->mom.hmax = 0;
}

References NODE::center, NODE::ChildType, NODE::cofm, NODE::DependsOnLocalMass, NODE::f, NODE::father, NODE::hmax, NODE::InternalTopLevel, NODE::len, NODE::mass, NODE::mom, NMAXCHILD, NodeChild::noccupied, PARTICLE_NODE_TYPE, NODE::s, NODE::sibling, NodeChild::suns, NODE::TopLevel, NodeChild::Types, and NODE::unused.

Referenced by create_new_node_layer(), and force_create_node_for_topnode().

Here is the caller graph for this function:

inside_node()

static int inside_node ( const struct NODE *  node, const int  p_i  )

inlinestatic

Definition at line 291 of file forcetree.c.

{
    /*One can also use a loop, but the compiler unrolls it only at -O3,
     *so this is a little faster*/
    int inside =
        (fabs(2*(P[p_i].Pos[0] - node->center[0])) <= node->len) *
        (fabs(2*(P[p_i].Pos[1] - node->center[1])) <= node->len) *
        (fabs(2*(P[p_i].Pos[2] - node->center[2])) <= node->len);
    return inside;
}

References NODE::center, NODE::len, and P.

Referenced by force_tree_create_nodes().

Here is the caller graph for this function:

merge_partial_force_trees()

int merge_partial_force_trees	(	int	left,
		int	right,
		struct NodeCache *	nc,
		int *	nnext,
		const struct ForceTree	tb,
		int	HybridNuGrav
	)

Definition at line 533 of file forcetree.c.

{
    int this_left = left;
    int this_right = right;
    const int left_end = tb.Nodes[left].sibling;
    const int right_end = tb.Nodes[right].sibling;
//     message(5, "Ends: %d %d\n", left_end, right_end);
    while(this_left != left_end && this_right != right_end)
    {
        if(this_left < tb.firstnode || this_right < tb.firstnode)
            endrun(10, "Encountered invalid node: %d %d < first %d\n", this_left, this_right, tb.firstnode);
        struct NODE * nleft = &tb.Nodes[this_left];
        struct NODE * nright = &tb.Nodes[this_right];
        if(nc->nnext_thread >= tb.lastnode)
            return 1;
#ifdef DEBUG
        /* Stop when we reach another topnode*/
        if((nleft->f.TopLevel && this_left != left) || (nright->f.TopLevel && this_right != right))
            endrun(6, "Encountered another topnode: left %d == right %d! type %d\n", this_left, this_right, nleft->f.ChildType);
        if(this_left == this_right)
            endrun(6, "Odd: left %d == right %d! type %d\n", this_left, this_right, nleft->f.ChildType);
//         message(1, "left %d right %d\n", this_left, this_right);
        /* Trees should be synced*/
        if(fabs(nleft->len / nright->len-1) > 1e-6)
            endrun(6, "Merge unsynced trees: %d %d len %g %g\n", this_left, this_right, nleft->len, nright->len);
#endif
        /* Two node nodes: keep walking down*/
        if(nleft->f.ChildType == NODE_NODE_TYPE && nright->f.ChildType == NODE_NODE_TYPE) {
            if(tb.Nodes[nleft->s.suns[0]].father < 0 || tb.Nodes[nright->s.suns[0]].father < 0)
                endrun(7, "Walking to nodes (%d %d) from (%d %d) fathers (%d %d)\n",
                       nleft->s.suns[0], nright->s.suns[0], this_left, this_right, tb.Nodes[nleft->s.suns[0]].father, tb.Nodes[nright->s.suns[0]].father);
            this_left = nleft->s.suns[0];
            this_right = nright->s.suns[0];
            continue;
        }
        /* If the right node has particles, add them to the left node, go to sibling on right and left.*/
        else if(nright->f.ChildType == PARTICLE_NODE_TYPE) {
            int i;
            for(i = 0; i < nright->s.noccupied; i++) {
                if(nright->s.suns[i] >= tb.firstnode)
                    endrun(8, "Bad child %d of %d\n", i, nright->s.suns[i], this_right);
                if(add_particle_to_tree(nright->s.suns[i], this_left, tb, HybridNuGrav, nc, nnext) < 0)
                    return 1;
            }
            /* Make sure that nodes which have
             * this_right as a sibling (there will
             * be a max of one, as it is a particle node
             * with no children) point to the replacement
             * on the left*/
            /* This condition is checking for the root node, which has no siblings*/
            if(this_right > right) {
                /* Find the father, then the next child*/
                struct NODE * fat = &tb.Nodes[nright->father];
                /* Find the position of this child in the father*/
                int sunloc = 0;
                for(i = 0; i < 8; i++)
                {
                    if(fat->s.suns[i] == this_right) {
                        sunloc = i;
                        break;
                    }
                }
                /* Change the sibling of the child next to this one*/
                if(sunloc > 0) {
                    if(tb.Nodes[fat->s.suns[i-1]].sibling == this_right)
                        tb.Nodes[fat->s.suns[i-1]].sibling = this_left;
                }
            }
            /* Mark the right node as now invalid*/
            nright->father = -5;
            /* Now go to sibling*/
            this_left = nleft->sibling;
            this_right = nright->sibling;
            continue;
        }
        /* If the left node has particles, add them to the right node,
         * then copy the right node over the left node and go to (old) sibling on right and left.*/
        else if(nleft->f.ChildType == PARTICLE_NODE_TYPE && nright->f.ChildType == NODE_NODE_TYPE) {
            /* Add the left particles to the right*/
            int i;
            for(i = 0; i < nleft->s.noccupied; i++) {
                if(nleft->s.suns[i] >= tb.firstnode)
                    endrun(8, "Bad child %d of %d\n", i, nleft->s.suns[i], this_left);
                if(add_particle_to_tree(nleft->s.suns[i], this_right, tb, HybridNuGrav, nc, nnext) < 0)
                    return 1;
            }
            /* Copy the right node over the left*/
            memmove(&nleft->s, &nright->s, sizeof(nleft->s));
            nleft->f.ChildType = NODE_NODE_TYPE;
            /* Zero the momenta for the parent*/
            memset(&nleft->mom, 0, sizeof(nleft->mom));
            /* Reset children to the new parent:
             * this assumes nright is a NODE NODE*/
            for(i = 0; i < 8; i++) {
                int child = nleft->s.suns[i];
                tb.Nodes[child].father = this_left;
            }
            /* Make sure final child points to the parent's sibling.*/
#ifdef DEBUG
            int oldsib = tb.Nodes[nleft->s.suns[7]].sibling;
#endif
            /* Walk downwards making sure all the children point to the new sibling.
             * Note also changes last particle node child. */
            int nn = this_left;
            while(tb.Nodes[nn].f.ChildType == NODE_NODE_TYPE) {
                nn = tb.Nodes[nn].s.suns[7];
#ifdef DEBUG
                if(tb.Nodes[nn].sibling != oldsib)
                    endrun(20, "Not the expected sibling %d != %d\n",tb.Nodes[nn].sibling, oldsib);
#endif
                tb.Nodes[nn].sibling = nleft->sibling;
            }
            /* Mark the right node as now invalid*/
            nright->father = -5;
            /* Next iteration is going to sibling*/
            this_left = nleft->sibling;
            this_right = nright->sibling;
            continue;
        }
        else
            endrun(6, "Nodes %d %d have unexpected type %d %d\n", this_left, this_right, nleft->f.ChildType, nright->f.ChildType);
    }
    return 0;
}

References add_particle_to_tree(), NODE::ChildType, endrun(), NODE::f, NODE::father, ForceTree::firstnode, ForceTree::lastnode, NODE::len, NODE::mom, NodeCache::nnext_thread, NodeChild::noccupied, NODE_NODE_TYPE, ForceTree::Nodes, PARTICLE_NODE_TYPE, NODE::s, NODE::sibling, NodeChild::suns, and NODE::TopLevel.

Referenced by force_tree_create_nodes().

Here is the call graph for this function:

Here is the caller graph for this function:

modify_internal_node()

static int modify_internal_node ( int  parent, int  subnode, int  p_toplace, const ForceTree  tb, const int  HybridNuGrav  )

static

Definition at line 361 of file forcetree.c.

{
    tb.Father[p_toplace] = parent;
    tb.Nodes[parent].s.suns[subnode] = p_toplace;
    /* Encode the type in the Types array*/
    tb.Nodes[parent].s.Types += P[p_toplace].Type << (3*subnode);
    if(!HybridNuGrav || P[p_toplace].Type != ForceTreeParams.FastParticleType)
        add_particle_moment_to_node(&tb.Nodes[parent], p_toplace);
    return 0;
}

References add_particle_moment_to_node(), forcetree_params::FastParticleType, ForceTree::Father, ForceTreeParams, ForceTree::Nodes, P, NODE::s, NodeChild::suns, and NodeChild::Types.

Referenced by add_particle_to_tree(), and create_new_node_layer().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

ForceTreeParams

struct forcetree_params ForceTreeParams

static

Referenced by force_tree_build(), init_forcetree_params(), and modify_internal_node().