treewalk.h File Reference

#include <stdint.h>
#include "utils/paramset.h"
#include "forcetree.h"

Include dependency graph for treewalk.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	TreeWalkQueryBase
struct	TreeWalkResultBase
struct	TreeWalkNgbIterBase
struct	LocalTreeWalk
struct	TreeWalk

Macros
#define	NODELISTLENGTH   8
#define	TREEWALK_REDUCE(A, B)   (A) = (mode==TREEWALK_PRIMARY)?(B):((A) + (B))
#define	MAXITER   400

Typedefs
typedef struct TreeWalk	TreeWalk
typedef int(*	TreeWalkVisitFunction) (TreeWalkQueryBase *input, TreeWalkResultBase *output, LocalTreeWalk *lv)
typedef void(*	TreeWalkNgbIterFunction) (TreeWalkQueryBase *input, TreeWalkResultBase *output, TreeWalkNgbIterBase *iter, LocalTreeWalk *lv)
typedef int(*	TreeWalkHasWorkFunction) (const int i, TreeWalk *tw)
typedef void(*	TreeWalkProcessFunction) (const int i, TreeWalk *tw)
typedef void(*	TreeWalkFillQueryFunction) (const int j, TreeWalkQueryBase *query, TreeWalk *tw)
typedef void(*	TreeWalkReduceResultFunction) (const int j, TreeWalkResultBase *result, const enum TreeWalkReduceMode mode, TreeWalk *tw)

Enumerations
enum	NgbTreeFindSymmetric { NGB_TREEFIND_SYMMETRIC , NGB_TREEFIND_ASYMMETRIC }
enum	TreeWalkReduceMode { TREEWALK_PRIMARY , TREEWALK_GHOSTS }
enum	TreeWalkType { TREEWALK_ACTIVE = 0 , TREEWALK_ALL , TREEWALK_SPLIT }

Functions
void	set_treewalk_params (ParameterSet *ps)
void	treewalk_run (TreeWalk *tw, int *active_set, size_t size)
int	treewalk_visit_ngbiter (TreeWalkQueryBase *I, TreeWalkResultBase *O, LocalTreeWalk *lv)
int	treewalk_export_particle (LocalTreeWalk *lv, int no)
int	treewalk_visit_nolist_ngbiter (TreeWalkQueryBase *I, TreeWalkResultBase *O, LocalTreeWalk *lv)
void	treewalk_do_hsml_loop (TreeWalk *tw, int *queue, int64_t queuesize, int update_hsml)
double	ngb_narrow_down (double *right, double *left, const double *radius, const double *numNgb, int maxcmpt, int desnumngb, int *closeidx, double BoxSize)
void	treewalk_build_queue (TreeWalk *tw, int *active_set, const size_t size, int may_have_garbage)

Macro Definition Documentation

MAXITER

#define MAXITER   400

Definition at line 201 of file treewalk.h.

NODELISTLENGTH

#define NODELISTLENGTH   8

Definition at line 8 of file treewalk.h.

TREEWALK_REDUCE

#define TREEWALK_REDUCE	(		A,
			B
	)		(A) = (mode==TREEWALK_PRIMARY)?(B):((A) + (B))

Definition at line 189 of file treewalk.h.

Typedef Documentation

TreeWalk

typedef struct TreeWalk TreeWalk

Definition at line 1 of file treewalk.h.

TreeWalkFillQueryFunction

typedef void(* TreeWalkFillQueryFunction) (const int j, TreeWalkQueryBase *query, TreeWalk *tw)

Definition at line 75 of file treewalk.h.

TreeWalkHasWorkFunction

typedef int(* TreeWalkHasWorkFunction) (const int i, TreeWalk *tw)

Definition at line 72 of file treewalk.h.

TreeWalkNgbIterFunction

typedef void(* TreeWalkNgbIterFunction) (TreeWalkQueryBase *input, TreeWalkResultBase *output, TreeWalkNgbIterBase *iter, LocalTreeWalk *lv)

Definition at line 70 of file treewalk.h.

TreeWalkProcessFunction

typedef void(* TreeWalkProcessFunction) (const int i, TreeWalk *tw)

Definition at line 73 of file treewalk.h.

TreeWalkReduceResultFunction

typedef void(* TreeWalkReduceResultFunction) (const int j, TreeWalkResultBase *result, const enum TreeWalkReduceMode mode, TreeWalk *tw)

Definition at line 76 of file treewalk.h.

TreeWalkVisitFunction

typedef int(* TreeWalkVisitFunction) (TreeWalkQueryBase *input, TreeWalkResultBase *output, LocalTreeWalk *lv)

Definition at line 68 of file treewalk.h.

Enumeration Type Documentation

NgbTreeFindSymmetric

enum NgbTreeFindSymmetric

Enumerator
NGB_TREEFIND_SYMMETRIC
NGB_TREEFIND_ASYMMETRIC

Definition at line 10 of file treewalk.h.

                          {
    NGB_TREEFIND_SYMMETRIC,
    NGB_TREEFIND_ASYMMETRIC,
};

TreeWalkReduceMode

enum TreeWalkReduceMode

Enumerator
TREEWALK_PRIMARY
TREEWALK_GHOSTS

Definition at line 15 of file treewalk.h.

                        {
    TREEWALK_PRIMARY,
    TREEWALK_GHOSTS,
};

TreeWalkType

enum TreeWalkType

Enumerator
TREEWALK_ACTIVE
TREEWALK_ALL
TREEWALK_SPLIT

Definition at line 78 of file treewalk.h.

                  {
    TREEWALK_ACTIVE = 0,
    TREEWALK_ALL,
    TREEWALK_SPLIT,
};

Function Documentation

ngb_narrow_down()

double ngb_narrow_down	(	double *	right,
		double *	left,
		const double *	radius,
		const double *	numNgb,
		int	maxcmpt,
		int	desnumngb,
		int *	closeidx,
		double	BoxSize
	)

Definition at line 1366 of file treewalk.c.

{
    int j;
    int close = 0;
    double ngbdist = fabs(numNgb[0] - desnumngb);
    for(j = 1; j < maxcmpt; j++){
        double newdist = fabs(numNgb[j] - desnumngb);
        if(newdist < ngbdist){
            ngbdist = newdist;
            close = j;
        }
    }
    if(closeidx)
        *closeidx = close;
 
    for(j = 0; j < maxcmpt; j++){
        if(numNgb[j] < desnumngb)
            *left = radius[j];
        if(numNgb[j] > desnumngb){
            *right = radius[j];
            break;
        }
    }
 
    double hsml = radius[close];
 
    if(*right > 0.99 * BoxSize){
        double dngbdv = 0;
        if(maxcmpt > 1 && (radius[maxcmpt-1]>radius[maxcmpt-2]))
            dngbdv = (numNgb[maxcmpt-1]-numNgb[maxcmpt-2])/(pow(radius[maxcmpt-1],3) - pow(radius[maxcmpt-2],3));
        /* Increase hsml by a maximum factor to avoid madness. We can be fairly aggressive about this factor.*/
        double newhsml = 4 * hsml;
        if(dngbdv > 0) {
            double dngb = (desnumngb - numNgb[maxcmpt-1]);
            double newvolume = pow(hsml,3) + dngb / dngbdv;
            if(pow(newvolume, 1./3) < newhsml)
                newhsml = pow(newvolume, 1./3);
        }
        hsml = newhsml;
    }
    if(hsml > *right)
        hsml = *right;
 
    if(*left == 0) {
        /* Extrapolate using volume, ie locally constant density*/
        double dngbdv = 0;
        if(radius[1] > radius[0])
            dngbdv = (numNgb[1] - numNgb[0]) / (pow(radius[1],3) - pow(radius[0],3));
        /* Derivative is not defined for minimum, so use 0.*/
        if(maxcmpt == 1 && radius[0] > 0)
            dngbdv = numNgb[0] / pow(radius[0],3);
 
        if(dngbdv > 0) {
            double dngb = desnumngb - numNgb[0];
            double newvolume = pow(hsml,3) + dngb / dngbdv;
            hsml = pow(newvolume, 1./3);
        }
    }
    if(hsml < *left)
        hsml = *left;
 
    return hsml;
}

Referenced by sfr_wind_weight_postprocess(), and stellar_density_check_neighbours().

Here is the caller graph for this function:

set_treewalk_params()

void set_treewalk_params

(

ParameterSet *

ps

)

Definition at line 55 of file treewalk.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0)
        ImportBufferBoost = param_get_int(ps, "ImportBufferBoost");
    MPI_Bcast(&ImportBufferBoost, 1, MPI_INT, 0, MPI_COMM_WORLD);
}

References ImportBufferBoost, param_get_int(), and ThisTask.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

treewalk_build_queue()

void treewalk_build_queue	(	TreeWalk *	tw,
		int *	active_set,
		const size_t	size,
		int	may_have_garbage
	)

Definition at line 394 of file treewalk.c.

{
    tw->NThread = omp_get_max_threads();
 
    if(!tw->haswork && !may_have_garbage)
    {
        tw->WorkSetSize = size;
        tw->WorkSet = active_set;
        tw->work_set_stolen_from_active = 1;
        return;
    }
 
    /* Since we use a static schedule below we only need size / tw->NThread elements per thread.
     * Add NThread so that each thread has a little extra space.*/
    size_t tsize = size / tw->NThread + tw->NThread;
    /*Watch out: tw->WorkSet may change a few lines later due to the realloc*/
    tw->WorkSet = (int *) mymalloc("ActiveQueue", tsize * sizeof(int) * tw->NThread);
    tw->work_set_stolen_from_active = 0;
    size_t i;
    size_t nqueue = 0;
 
    /*We want a lockless algorithm which preserves the ordering of the particle list.*/
    size_t *nqthr = ta_malloc("nqthr", size_t, tw->NThread);
    int **thrqueue = ta_malloc("thrqueue", int *, tw->NThread);
 
    gadget_setup_thread_arrays(tw->WorkSet, thrqueue, nqthr, tsize, tw->NThread);
 
    /* We enforce schedule static to ensure that each thread executes on contiguous particles.
     * Note static enforces the monotonic modifier but on OpenMP 5.0 nonmonotonic is the default.
     * static also ensures that no single thread gets more than tsize elements.*/
    size_t schedsz = size/tw->NThread+1;
    #pragma omp parallel for schedule(static, schedsz)
    for(i=0; i < size; i++)
    {
        const int tid = omp_get_thread_num();
        /*Use raw particle number if active_set is null, otherwise use active_set*/
        const int p_i = active_set ? active_set[i] : (int) i;
 
        /* Skip the garbage particles */
        if(P[p_i].IsGarbage)
            continue;
 
        if(tw->haswork && !tw->haswork(p_i, tw))
            continue;
#ifdef DEBUG
        if(nqthr[tid] >= tsize)
            endrun(5, "tid = %d nqthr = %d, tsize = %d size = %d, tw->Nthread = %d i = %d\n", tid, nqthr[tid], tsize, size, tw->NThread, i);
#endif
        thrqueue[tid][nqthr[tid]] = p_i;
        nqthr[tid]++;
    }
    /*Merge step for the queue.*/
    nqueue = gadget_compact_thread_arrays(tw->WorkSet, thrqueue, nqthr, tw->NThread);
    ta_free(thrqueue);
    ta_free(nqthr);
    /*Shrink memory*/
    tw->WorkSet = (int *) myrealloc(tw->WorkSet, sizeof(int) * nqueue);
 
#if 0
    /* check the uniqueness of the active_set list. This is very slow. */
    qsort_openmp(tw->WorkSet, nqueue, sizeof(int), cmpint);
    for(i = 0; i < nqueue - 1; i ++) {
        if(tw->WorkSet[i] == tw->WorkSet[i+1]) {
            endrun(8829, "A few particles are twicely active.\n");
        }
    }
#endif
    tw->WorkSetSize = nqueue;
}

References endrun(), gadget_compact_thread_arrays(), gadget_setup_thread_arrays(), TreeWalk::haswork, mymalloc, myrealloc, TreeWalk::NThread, P, qsort_openmp, ta_free, ta_malloc, TreeWalk::work_set_stolen_from_active, TreeWalk::WorkSet, and TreeWalk::WorkSetSize.

Referenced by blackhole(), ev_begin(), and treewalk_do_hsml_loop().

Here is the call graph for this function:

Here is the caller graph for this function:

treewalk_do_hsml_loop()

void treewalk_do_hsml_loop	(	TreeWalk *	tw,
		int *	queue,
		int64_t	queuesize,
		int	update_hsml
	)

Definition at line 1254 of file treewalk.c.

{
    int64_t ntot = 0;
    int NumThreads = omp_get_max_threads();
    tw->NPLeft = ta_malloc("NPLeft", size_t, NumThreads);
    tw->NPRedo = ta_malloc("NPRedo", int *, NumThreads);
    tw->maxnumngb = ta_malloc("numngb", double, NumThreads);
    tw->minnumngb = ta_malloc("numngb2", double, NumThreads);
 
    /* Build the first queue */
    treewalk_build_queue(tw, queue, queuesize, 0);
    /* Next call to treewalk_run will over-write these pointers*/
    int64_t size = tw->WorkSetSize;
    int * ReDoQueue = tw->WorkSet;
    /* First queue is allocated low*/
    int alloc_high = 0;
    /* We don't need to redo the queue generation
     * but need to keep track of allocated memory.*/
    int orig_queue_alloc = (tw->haswork != NULL);
    tw->haswork = NULL;
 
    /* we will repeat the whole thing for those particles where we didn't find enough neighbours */
    do {
        /* The RedoQueue needs enough memory to store every workset particle on every thread, because
         * we cannot guarantee that the sph particles are evenly spread across threads!*/
        int * CurQueue = ReDoQueue;
        int i;
        for(i = 0; i < NumThreads; i++) {
            tw->maxnumngb[i] = 0;
            tw->minnumngb[i] = 1e50;
        }
 
        /* The ReDoQueue swaps between high and low allocations so we can have two allocated alternately*/
        if(update_hsml) {
            if(!alloc_high) {
                ReDoQueue = (int *) mymalloc2("ReDoQueue", size * sizeof(int) * NumThreads);
                alloc_high = 1;
            }
            else {
                ReDoQueue = (int *) mymalloc("ReDoQueue", size * sizeof(int) * NumThreads);
                alloc_high = 0;
            }
            tw->Redo_thread_alloc = size;
            gadget_setup_thread_arrays(ReDoQueue, tw->NPRedo, tw->NPLeft, size, NumThreads);
        }
        treewalk_run(tw, CurQueue, size);
 
        /* Now done with the current queue*/
        if(orig_queue_alloc || tw->Niteration > 1)
            myfree(CurQueue);
 
        /* We can stop if we are not updating hsml*/
        if(!update_hsml)
            break;
 
        /* Set up the next queue*/
        size = gadget_compact_thread_arrays(ReDoQueue, tw->NPRedo, tw->NPLeft, NumThreads);
 
        MPI_Allreduce(&size, &ntot, 1, MPI_INT64, MPI_SUM, MPI_COMM_WORLD);
        if(ntot == 0){
            myfree(ReDoQueue);
            break;
        }
        for(i = 1; i < NumThreads; i++) {
            if(tw->maxnumngb[0] < tw->maxnumngb[i])
                tw->maxnumngb[0] = tw->maxnumngb[i];
            if(tw->minnumngb[0] > tw->minnumngb[i])
                tw->minnumngb[0] = tw->minnumngb[i];
        }
        double minngb, maxngb;
        MPI_Reduce(&tw->maxnumngb[0], &maxngb, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
        MPI_Reduce(&tw->minnumngb[0], &minngb, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
        message(0, "Max ngb=%g, min ngb=%g\n", maxngb, minngb);
 
        /*Shrink memory*/
        ReDoQueue = (int *) myrealloc(ReDoQueue, sizeof(int) * size);
 
        /*
        if(ntot < 1 ) {
            foreach(ActiveParticle)
            {
                if(density_haswork(i)) {
                    MyFloat Left = DENSITY_GET_PRIV(tw)->Left[i];
                    MyFloat Right = DENSITY_GET_PRIV(tw)->Right[i];
                    message (1, "i=%d task=%d ID=%llu type=%d, Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g\n   pos=(%g|%g|%g)\n",
                         i, ThisTask, P[i].ID, P[i].Type, P[i].Hsml, Left, Right,
                         (float) P[i].NumNgb, Right - Left, P[i].Pos[0], P[i].Pos[1], P[i].Pos[2]);
                }
            }
 
        }
        */
#ifdef DEBUG
        if(ntot == 1 && size > 0 && tw->Niteration > 20 ) {
            int pp = ReDoQueue[0];
            message(1, "Remaining i=%d, t %d, pos %g %g %g, hsml: %g\n", pp, P[pp].Type, P[pp].Pos[0], P[pp].Pos[1], P[pp].Pos[2], P[pp].Hsml);
        }
#endif
 
        if(size > 0 && tw->Niteration > MAXITER) {
            endrun(1155, "failed to converge density for %ld particles\n", ntot);
        }
    } while(1);
    ta_free(tw->minnumngb);
    ta_free(tw->maxnumngb);
    ta_free(tw->NPRedo);
    ta_free(tw->NPLeft);
}

References endrun(), gadget_compact_thread_arrays(), gadget_setup_thread_arrays(), TreeWalk::haswork, MAXITER, TreeWalk::maxnumngb, message(), TreeWalk::minnumngb, MPI_INT64, myfree, mymalloc, mymalloc2, myrealloc, TreeWalk::Niteration, TreeWalk::NPLeft, TreeWalk::NPRedo, P, TreeWalk::Redo_thread_alloc, ta_free, ta_malloc, treewalk_build_queue(), treewalk_run(), TreeWalk::WorkSet, and TreeWalk::WorkSetSize.

Referenced by density(), stellar_density(), and winds_find_weights().

Here is the call graph for this function:

Here is the caller graph for this function:

treewalk_export_particle()

int treewalk_export_particle	(	LocalTreeWalk *	lv,
		int	no
	)

Definition at line 567 of file treewalk.c.

{
    if(lv->mode != 0) {
        endrun(1, "Trying to export a ghost particle.\n");
    }
    const int target = lv->target;
    int *exportflag = lv->exportflag;
    int *exportnodecount = lv->exportnodecount;
    size_t *exportindex = lv->exportindex;
    TreeWalk * tw = lv->tw;
 
    const int task = tw->tree->TopLeaves[no - tw->tree->lastnode].Task;
 
    if(exportflag[task] != target)
    {
        exportflag[task] = target;
        exportnodecount[task] = NODELISTLENGTH;
    }
 
    if(exportnodecount[task] == NODELISTLENGTH)
    {
        /* out of buffer space. Need to interrupt. */
        if(lv->Nexport >= lv->BunchSize) {
            return -1;
        }
        /* This index is a unique entry in the global DataNodeList and DataIndexTable.*/
        size_t nexp = lv->Nexport + lv->DataIndexOffset;
        exportnodecount[task] = 0;
        exportindex[task] = nexp;
        DataIndexTable[nexp].Task = task;
        DataIndexTable[nexp].Index = target;
        DataIndexTable[nexp].IndexGet = nexp;
        lv->Nexport++;
        lv->NThisParticleExport++;
    }
 
    /* Set the NodeList entry*/
    DataNodeList[exportindex[task]].NodeList[exportnodecount[task]++] =
            tw->tree->TopLeaves[no - tw->tree->lastnode].treenode;
 
    if(exportnodecount[task] < NODELISTLENGTH)
        DataNodeList[exportindex[task]].NodeList[exportnodecount[task]] = -1;
    return 0;
}

References LocalTreeWalk::BunchSize, LocalTreeWalk::DataIndexOffset, DataIndexTable, DataNodeList, endrun(), LocalTreeWalk::exportflag, LocalTreeWalk::exportindex, LocalTreeWalk::exportnodecount, data_index::Index, data_index::IndexGet, ForceTree::lastnode, LocalTreeWalk::mode, LocalTreeWalk::Nexport, data_nodelist::NodeList, NODELISTLENGTH, LocalTreeWalk::NThisParticleExport, LocalTreeWalk::target, topleaf_data::Task, data_index::Task, ForceTree::TopLeaves, TreeWalk::tree, topleaf_data::treenode, and LocalTreeWalk::tw.

Referenced by force_treeev_shortrange(), ngb_treefind_threads(), and treewalk_visit_nolist_ngbiter().

Here is the call graph for this function:

Here is the caller graph for this function:

treewalk_run()

void treewalk_run	(	TreeWalk *	tw,
		int *	active_set,
		size_t	size
	)

Definition at line 619 of file treewalk.c.

{
    if(!force_tree_allocated(tw->tree)) {
        endrun(0, "Tree has been freed before this treewalk.\n");
    }
 
    GDB_current_ev = tw;
 
    ev_begin(tw, active_set, size);
 
    if(tw->preprocess) {
        int64_t i;
        #pragma omp parallel for
        for(i = 0; i < tw->WorkSetSize; i ++) {
            const int p_i = tw->WorkSet ? tw->WorkSet[i] : i;
            tw->preprocess(p_i, tw);
        }
    }
 
    if(tw->visit) {
        tw->Nexportfull = 0;
        tw->evaluated = NULL;
        do
        {
            /* Keep track of which particles have been evaluated across buffer fill ups.
             * Do this if we are not allowed to evaluate anything twice,
             * or if the buffer filled up already.*/
            if((!tw->evaluated) && (tw->Nexportfull == 1 || tw->repeatdisallowed)) {
                tw->evaluated = (char *) mymalloc("evaluated", sizeof(char)*tw->WorkSetSize);
                memset(tw->evaluated, 0, sizeof(char)*tw->WorkSetSize);
            }
            ev_primary(tw); /* do local particles and prepare export list */
            /* exchange particle data */
            const struct SendRecvBuffer sndrcv = ev_get_remote(tw);
            /* now do the particles that were sent to us */
            ev_secondary(tw);
 
            /* import the result to local particles */
            ev_reduce_result(sndrcv, tw);
 
            tw->Nexportfull ++;
            tw->Nexport_sum += tw->Nexport;
            ta_free(sndrcv.Send_count);
        } while(ev_ndone(tw) < tw->NTask);
        if(tw->evaluated)
            myfree(tw->evaluated);
    }
 
#ifdef DEBUG
    /*int64_t totNodesinlist, totlist;
    MPI_Reduce(&tw->Nnodesinlist,  &totNodesinlist, 1, MPI_INT64, MPI_SUM, 0, MPI_COMM_WORLD);
    MPI_Reduce(&tw->Nlist,  &totlist, 1, MPI_INT64, MPI_SUM, 0, MPI_COMM_WORLD);
    message(0, "Nodes in nodelist: %g (avg). %ld nodes, %ld lists\n", ((double) totNodesinlist)/totlist, totlist, totNodesinlist);*/
#endif
    double tstart, tend;
 
    tstart = second();
 
    if(tw->postprocess) {
        int64_t i;
        #pragma omp parallel for
        for(i = 0; i < tw->WorkSetSize; i ++) {
            const int p_i = tw->WorkSet ? tw->WorkSet[i] : i;
            tw->postprocess(p_i, tw);
        }
    }
    tend = second();
    tw->timecomp3 = timediff(tstart, tend);
    ev_finish(tw);
    tw->Niteration++;
}

References endrun(), ev_begin(), ev_finish(), ev_get_remote(), ev_ndone(), ev_primary(), ev_reduce_result(), ev_secondary(), TreeWalk::evaluated, force_tree_allocated(), GDB_current_ev, myfree, mymalloc, TreeWalk::Nexport, TreeWalk::Nexport_sum, TreeWalk::Nexportfull, TreeWalk::Niteration, TreeWalk::NTask, TreeWalk::postprocess, TreeWalk::preprocess, TreeWalk::repeatdisallowed, second(), SendRecvBuffer::Send_count, ta_free, TreeWalk::timecomp3, timediff(), TreeWalk::tree, TreeWalk::visit, TreeWalk::WorkSet, and TreeWalk::WorkSetSize.

Referenced by blackhole(), fof_label_primary(), fof_label_secondary(), grav_short_pair(), grav_short_tree(), hydro_force(), ionize_all_part(), metal_return(), treewalk_do_hsml_loop(), and winds_and_feedback().

Here is the call graph for this function:

Here is the caller graph for this function:

treewalk_visit_ngbiter()

int treewalk_visit_ngbiter	(	TreeWalkQueryBase *	I,
		TreeWalkResultBase *	O,
		LocalTreeWalk *	lv
	)

Definition at line 924 of file treewalk.c.

{
 
    TreeWalkNgbIterBase * iter = (TreeWalkNgbIterBase *) alloca(lv->tw->ngbiter_type_elsize);
 
    /* Kick-start the iteration with other == -1 */
    iter->other = -1;
    lv->tw->ngbiter(I, O, iter, lv);
#ifdef DEBUG
    /* Check whether the tree contains the particles we are looking for*/
    if((lv->tw->tree->mask & iter->mask) == 0)
        endrun(5, "Tried to run treewalk for particle mask %d but tree mask is %d.\n", iter->mask, lv->tw->tree->mask);
#endif
    const double BoxSize = lv->tw->tree->BoxSize;
 
    int ninteractions = 0;
    int inode = 0;
 
    for(inode = 0; (lv->mode == 0 && inode < 1)|| (lv->mode == 1 && inode < NODELISTLENGTH && I->NodeList[inode] >= 0); inode++)
    {
        int numcand = ngb_treefind_threads(I, O, iter, I->NodeList[inode], lv);
        /* Export buffer is full end prematurally */
        if(numcand < 0) return numcand;
 
        /* If we are here, export is succesful. Work on the this particle -- first
         * filter out all of the candidates that are actually outside. */
        int numngb;
 
        for(numngb = 0; numngb < numcand; numngb ++) {
            int other = lv->ngblist[numngb];
 
            /* Skip garbage*/
            if(P[other].IsGarbage)
                continue;
            /* In case the type of the particle has changed since the tree was built.
             * Happens for wind treewalk for gas turned into stars on this timestep.*/
            if(!((1<<P[other].Type) & iter->mask)) {
                continue;
            }
 
            double dist;
 
            if(iter->symmetric == NGB_TREEFIND_SYMMETRIC) {
                dist = DMAX(P[other].Hsml, iter->Hsml);
            } else {
                dist = iter->Hsml;
            }
 
            double r2 = 0;
            int d;
            double h2 = dist * dist;
            for(d = 0; d < 3; d ++) {
                /* the distance vector points to 'other' */
                iter->dist[d] = NEAREST(I->Pos[d] - P[other].Pos[d], BoxSize);
                r2 += iter->dist[d] * iter->dist[d];
                if(r2 > h2) break;
            }
            if(r2 > h2) continue;
 
            /* update the iter and call the iteration function*/
            iter->r2 = r2;
            iter->r = sqrt(r2);
            iter->other = other;
 
            lv->tw->ngbiter(I, O, iter, lv);
        }
 
        ninteractions += numngb;
    }
 
    lv->Ninteractions += ninteractions;
    if(lv->mode == 1) {
        lv->Nnodesinlist += inode;
        lv->Nlist += 1;
    }
    return 0;
}

References ForceTree::BoxSize, TreeWalkNgbIterBase::dist, DMAX, endrun(), TreeWalkNgbIterBase::Hsml, ForceTree::mask, TreeWalkNgbIterBase::mask, LocalTreeWalk::mode, NEAREST, NGB_TREEFIND_SYMMETRIC, ngb_treefind_threads(), TreeWalk::ngbiter, TreeWalk::ngbiter_type_elsize, LocalTreeWalk::ngblist, LocalTreeWalk::Ninteractions, LocalTreeWalk::Nlist, LocalTreeWalk::Nnodesinlist, TreeWalkQueryBase::NodeList, TreeWalkNgbIterBase::other, P, TreeWalkQueryBase::Pos, TreeWalkNgbIterBase::r, TreeWalkNgbIterBase::r2, TreeWalkNgbIterBase::symmetric, TreeWalk::tree, and LocalTreeWalk::tw.

Referenced by blackhole(), fof_label_primary(), grav_short_pair(), hydro_force(), ionize_all_part(), and metal_return().

Here is the call graph for this function:

Here is the caller graph for this function:

treewalk_visit_nolist_ngbiter()

int treewalk_visit_nolist_ngbiter	(	TreeWalkQueryBase *	I,
		TreeWalkResultBase *	O,
		LocalTreeWalk *	lv
	)

Definition at line 1143 of file treewalk.c.

{
    TreeWalkNgbIterBase * iter = (TreeWalkNgbIterBase *) alloca(lv->tw->ngbiter_type_elsize);
 
    /* Kick-start the iteration with other == -1 */
    iter->other = -1;
    lv->tw->ngbiter(I, O, iter, lv);
 
    int inode;
    for(inode = 0; (lv->mode == 0 && inode < 1)|| (lv->mode == 1 && inode < NODELISTLENGTH && I->NodeList[inode] >= 0); inode++)
    {
        int no = I->NodeList[inode];
        const ForceTree * tree = lv->tw->tree;
        const double BoxSize = tree->BoxSize;
 
        while(no >= 0)
        {
            struct NODE *current = &tree->Nodes[no];
 
            /* When walking exported particles we start from the encompassing top-level node,
            * so if we get back to a top-level node again we are done.*/
            if(lv->mode == 1) {
                /* The first node is always top-level*/
                if(current->f.TopLevel && no != I->NodeList[inode]) {
                    /* we reached a top-level node again, which means that we are done with the branch */
                    break;
                }
            }
 
            /* Cull the node */
            if(0 == cull_node(I, iter, current, BoxSize)) {
                /* in case the node can be discarded */
                no = current->sibling;
                continue;
            }
 
            /* Node contains relevant particles, add them.*/
            if(current->f.ChildType == PARTICLE_NODE_TYPE) {
                int i;
                int * suns = current->s.suns;
                for (i = 0; i < current->s.noccupied; i++) {
                    /* must be the correct type: compare the
                    * current type for this subnode extracted
                    * from the bitfield to the mask.*/
                    int type = (current->s.Types >> (3*i)) % 8;
 
                    if(!((1<<type) & iter->mask))
                        continue;
 
                    /* Now evaluate a particle for the list*/
                    int other = suns[i];
                    /* Skip garbage*/
                    if(P[other].IsGarbage)
                        continue;
                    /* In case the type of the particle has changed since the tree was built.
                    * Happens for wind treewalk for gas turned into stars on this timestep.*/
                    if(!((1<<P[other].Type) & iter->mask))
                        continue;
 
                    double dist = iter->Hsml;
                    double r2 = 0;
                    int d;
                    double h2 = dist * dist;
                    for(d = 0; d < 3; d ++) {
                        /* the distance vector points to 'other' */
                        iter->dist[d] = NEAREST(I->Pos[d] - P[other].Pos[d], BoxSize);
                        r2 += iter->dist[d] * iter->dist[d];
                        if(r2 > h2) break;
                    }
                    if(r2 > h2) continue;
 
                    /* update the iter and call the iteration function*/
                    iter->r2 = r2;
                    iter->other = other;
                    iter->r = sqrt(r2);
                    lv->tw->ngbiter(I, O, iter, lv);
                }
                /* Move sideways*/
                no = current->sibling;
                continue;
            }
            else if(current->f.ChildType == PSEUDO_NODE_TYPE) {
                /* pseudo particle */
                if(lv->mode == 1) {
                    endrun(12312, "Secondary for particle %d from node %d found pseudo at %d.\n", lv->target, I->NodeList[inode], current);
                } else {
                    /* Export the pseudo particle*/
                    if(-1 == treewalk_export_particle(lv, current->s.suns[0]))
                        return -1;
                    /* Move sideways*/
                    no = current->sibling;
                    continue;
                }
            }
            /* ok, we need to open the node */
            no = current->s.suns[0];
        }
    }
 
    if(lv->mode == 1) {
        lv->Nnodesinlist += inode;
        lv->Nlist += 1;
    }
    return 0;
}

References ForceTree::BoxSize, NODE::ChildType, cull_node(), TreeWalkNgbIterBase::dist, endrun(), NODE::f, TreeWalkNgbIterBase::Hsml, TreeWalkNgbIterBase::mask, LocalTreeWalk::mode, NEAREST, TreeWalk::ngbiter, TreeWalk::ngbiter_type_elsize, LocalTreeWalk::Nlist, LocalTreeWalk::Nnodesinlist, NodeChild::noccupied, TreeWalkQueryBase::NodeList, ForceTree::Nodes, TreeWalkNgbIterBase::other, P, PARTICLE_NODE_TYPE, TreeWalkQueryBase::Pos, PSEUDO_NODE_TYPE, TreeWalkNgbIterBase::r, TreeWalkNgbIterBase::r2, NODE::s, NODE::sibling, NodeChild::suns, LocalTreeWalk::target, NODE::TopLevel, TreeWalk::tree, treewalk_export_particle(), LocalTreeWalk::tw, and NodeChild::Types.

Referenced by density(), fof_label_secondary(), stellar_density(), and winds_find_weights().

Here is the call graph for this function:

Here is the caller graph for this function: