#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
Go to the source code of this file.
Defines  
#define  NTAB 1000 
#define  NEAREST(x) (((x)>boxhalf)?((x)boxsize):(((x)<boxhalf)?((x)+boxsize):(x))) 
#define  EN 64 
Functions  
int  force_treebuild (int npart) 
int  force_treebuild_single (int npart) 
void  force_create_empty_nodes (int no, int topnode, int bits, int x, int y, int z, int *nodecount, int *nextfree) 
void  force_insert_pseudo_particles (void) 
void  force_update_node_recursive (int no, int sib, int father) 
void  force_update_pseudoparticles (void) 
void  force_exchange_pseudodata (void) 
void  force_treeupdate_pseudos (void) 
void  force_flag_localnodes (void) 
void  force_update_len (void) 
void  force_update_node_len_local (void) 
void  force_update_node_len_toptree (void) 
void  force_update_hmax (void) 
void  force_update_node_hmax_local (void) 
void  force_update_node_hmax_toptree (void) 
int  force_treeevaluate (int target, int mode, double *ewaldcountsum) 
int  force_treeevaluate_shortrange (int target, int mode) 
int  force_treeevaluate_ewald_correction (int target, int mode, double pos_x, double pos_y, double pos_z, double aold) 
void  force_treeevaluate_potential (int target, int mode) 
void  force_treeevaluate_potential_shortrange (int target, int mode) 
void  force_treeallocate (int maxnodes, int maxpart) 
void  force_treefree (void) 
int  force_treeevaluate_direct (int target, int mode) 
void  dump_particles (void) 
void  ewald_init (void) 
void  ewald_corr (double dx, double dy, double dz, double *fper) 
double  ewald_pot_corr (double dx, double dy, double dz) 
double  ewald_psi (double x[3]) 
void  ewald_force (int iii, int jjj, int kkk, double x[3], double force[3]) 
This file contains the computation of the gravitational force by means of a tree. The type of tree implemented is a geometrical octtree, starting from a cube encompassing all particles. This cube is automatically found in the domain decomposition, which also splits up the global "toplevel" tree along node boundaries, moving the particles of different parts of the tree to separate processors. Tree nodes can be dynamically updated in drift/kick operations to avoid having to reconstruct the tree every timestep.
Definition in file forcetree.c.

Size of 3D lockup table for Ewald correction force Definition at line 45 of file forcetree.c. Referenced by ewald_pot_corr(), and force_treeevaluate_ewald_correction(). 

Macro that maps a distance to the nearest periodic neighbour Definition at line 43 of file forcetree.c. Referenced by force_treeevaluate_ewald_correction(). 

length of lockup table for shortrange force kernel in TreePM algorithm Definition at line 31 of file forcetree.c. 

This function dumps some of the basic particle data to a file. In case the tree construction fails, it is called just before the run terminates with an error message. Examination of the generated file may then give clues to what caused the problem. Definition at line 2763 of file forcetree.c. Referenced by force_create_empty_nodes(), and force_treebuild_single(). 

This function looks up the correction force due to the infinite number of periodic particle/node images. We here use trilinear interpolation to get it from the precomputed tables, which contain one octant around the target particle at the origin. The other octants are obtained from it by exploiting the symmetry properties. Definition at line 2951 of file forcetree.c. 

This function computes the force correction term (difference between full force of infinite lattice and nearest image) by Ewald summation. Definition at line 3134 of file forcetree.c. 

This function initializes tables with the correction force and the correction potential due to the periodic images of a point mass located at the origin. These corrections are obtained by Ewald summation. (See Hernquist, Bouchet, Suto, ApJS, 1991, 75, 231) The correction fields are used to obtain the full periodic force if periodic boundaries combined with the pure tree algorithm are used. For the TreePM algorithm, the Ewald correction is not used. The correction fields are stored on disk once they are computed. If a corresponding file is found, they are loaded from disk to speed up the initialization. The Ewald summation is done in parallel, i.e. the processors share the work to compute the tables if needed. Definition at line 2802 of file forcetree.c. References ThisTask. Referenced by begrun(). 

This function looks up the correction potential due to the infinite number of periodic particle/node images. We here use trilinear interpolation to get it from the precomputed table, which contains one octant around the target particle at the origin. The other octants are obtained from it by exploiting symmetry properties. Definition at line 3040 of file forcetree.c. References EN. 

This function computes the potential correction term by means of Ewald summation. Definition at line 3093 of file forcetree.c. 

This function recursively creates a set of empty tree nodes which corresponds to the toplevel tree for the domain grid. This is done to ensure that this toplevel tree is always "complete" so that we can easily associate the pseudoparticles of other CPUs with treenodes 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 294 of file forcetree.c. References NODE::center, topnode_data::Daughter, DomainNodeIndex, dump_particles(), endrun(), topnode_data::Leaf, NODE::len, MaxNodes, Nodes, peano_hilbert_key(), NODE::suns, ThisTask, TopNodes, and NODE::u. Referenced by force_treebuild_single(). 

This function communicates the values of the multipole moments of the toplevel treenodes of the domain grid. This data can then be used to update the pseudoparticles on each CPU accordingly. Definition at line 687 of file forcetree.c. References NODE::bitflags, NODE::d, DomainMoment, DomainNodeIndex, Extnodes, DomainNODE::mass, NODE::mass, Nodes, DomainNODE::s, NODE::s, NODE::u, DomainNODE::vs, and extNODE::vs. 

This function flags nodes in the toplevel tree that are dependent on local particle data. Definition at line 819 of file forcetree.c. References NODE::bitflags, NODE::d, DomainNodeIndex, NODE::father, Nodes, and NODE::u. Referenced by force_treebuild(). 

this function inserts pseudoparticles which will represent the mass distribution of the other CPUs. Initially, the mass of the pseudoparticles is set to zero, and their coordinate is set to the center of the domaincell they correspond to. These quantities will be updated later on. Definition at line 348 of file forcetree.c. References NODE::center, DomainMoment, DomainNodeIndex, DomainNODE::mass, Nodes, and DomainNODE::s. Referenced by force_treebuild_single(). 

This function allocates the memory used for storage of the tree and of auxiliary arrays needed for treewalk and linklists. Usually, maxnodes approximately equal to 0.7*maxpart is sufficient to store the tree for up to maxpart particles. < Maximum number of nodes in the toplevel tree used for domain decomposition < Maximum number of nodes in the toplevel tree used for domain decomposition Definition at line 2551 of file forcetree.c. References endrun(), MaxNodes, and Nodes_base. Referenced by init(), pmforce_periodic(), pmpotential_periodic(), and restart(). 

This function is a driver routine for constructing the gravitational octtree, which is done by calling a small number of other functions. Definition at line 61 of file forcetree.c. References All, force_flag_localnodes(), force_treebuild_single(), force_update_pseudoparticles(), Numnodestree, global_data_all_processes::Time, and TimeOfLastTreeConstruction. Referenced by compute_potential(), gravity_tree(), and ngb_treebuild(). 

Constructs the gravitational octtree. The index convention for accessing tree nodes is the following: the indices 0...NumPart1 reference single particles, the indices All.MaxPart.... All.MaxPart+nodes1 reference tree nodes. `Nodes_base' points to the first tree node, while `nodes' is shifted such that nodes[All.MaxPart] gives the first tree node. Finally, node indices with values 'All.MaxPart + MaxNodes' and larger indicate "pseudo particles", i.e. multipole moments of toplevel nodes that lie on different CPUs. If such a node needs to be opened, the corresponding particle must be exported to that CPU. The 'Extnodes' structure parallels that of 'Nodes'. Its information is only needed for the SPH part of the computation. (The data is split onto these two structures as a tuning measure. If it is merged into 'Nodes' a somewhat bigger size of the nodes also for gravity would result, which would reduce cache utilization slightly. < Bits per dimension available for PeanoHilbert order. Note: If peanokey is defined as type int, the allowed maximum is 10. If 64bit integers are used, the maximum is 21 Definition at line 93 of file forcetree.c. References All, BITS_PER_DIMENSION, NODE::center, NODE::d, topnode_data::Daughter, DomainCenter, DomainCorner, DomainFac, DomainNodeIndex, dump_particles(), endrun(), force_create_empty_nodes(), force_insert_pseudo_particles(), force_update_node_recursive(), global_data_all_processes::ForceSoftening, get_random_number(), particle_data::GravCost, topnode_data::Leaf, NODE::len, MaxNodes, global_data_all_processes::MaxPart, Nextnode, NODE::nextnode, Nodes, P, peano_hilbert_key(), peanokey, particle_data::Pos, topnode_data::Size, topnode_data::StartKey, NODE::suns, ThisTask, TopNodes, particle_data::Type, and NODE::u. Referenced by force_treebuild(). 

This routine computes the gravitational force for a given local particle, or for a particle in the communication buffer. Depending on the value of TypeOfOpeningCriterion, either the geometrical BH cellopening criterion, or the `relative' opening criterion is used. Definition at line 1111 of file forcetree.c. References All, global_data_all_processes::ErrTolForceAcc, sph_particle_data::Hsml, particle_data::OldAcc, P, particle_data::Pos, SphP, and particle_data::Type. Referenced by gravity_tree(). 

This function does the force computation with direct summation for the specified particle in the communication buffer. This can be useful for debugging purposes, in particular for explicit checks of the force accuracy. Definition at line 2632 of file forcetree.c. References All, P, particle_data::Pos, and particle_data::Type. Referenced by gravity_forcetest(). 

This function computes the Ewald correction, and is needed if periodic boundary conditions together with a pure tree algorithm are used. Note that the ordinary tree walk does not carry out this correction directly as it was done in Gadget1.1. Instead, the tree is walked a second time. This is actually faster because the "EwaldTreewalk" can use a different opening criterion than the normal tree walk. In particular, the Ewald correction is negligible for particles that are very close, but it is large for particles that are far away (this is quite different for the normal direct force). So we can here use a different opening criterion. Sufficient accuracy is usually obtained if the node length has dropped to a certain fraction ~< 0.25 of the BoxLength. However, we may only shortcut the interaction list of the normal full Ewald tree walk if we are sure that the whole node and all daughter nodes "lie on the same side" of the periodic boundary, i.e. that the real tree walk would not find a daughter node or particle that was mapped to a different nearest neighbour position when the tree walk would be further refined. Definition at line 1734 of file forcetree.c. References gravdata_in::Acc, All, NODE::bitflags, global_data_all_processes::BoxSize, NODE::center, NODE::d, DomainTask, EN, global_data_all_processes::ErrTolTheta, Exportflag, particle_data::GravAccel, particle_data::GravCost, GravDataResult, NODE::len, particle_data::Mass, NODE::mass, global_data_all_processes::MaxPart, NEAREST, Nextnode, NODE::nextnode, gravdata_in::Ninteractions, Nodes, P, particle_data::Pos, NODE::s, NODE::sibling, NODE::u, gravdata_in::u, and gravdata_in::w. 

This routine computes the gravitational potential by walking the tree. The same opening criteria is used as for the gravitational force walk. Definition at line 2000 of file forcetree.c. References All, global_data_all_processes::ErrTolForceAcc, sph_particle_data::Hsml, particle_data::OldAcc, P, particle_data::Pos, SphP, and particle_data::Type. Referenced by compute_potential(). 

This function computes the shortrange potential when the TreePM algorithm is used. This potential is the Newtonian potential, modified by a complementary error function. Definition at line 2247 of file forcetree.c. References All, global_data_all_processes::ErrTolForceAcc, sph_particle_data::Hsml, particle_data::OldAcc, P, particle_data::Pos, SphP, and particle_data::Type. Referenced by compute_potential(). 

In the TreePM algorithm, the tree is walked only locally around the target coordinate. Tree nodes that fall outside a box of half sidelength Rcut= RCUT*ASMTH*MeshSize can be discarded. The shortrange potential is modified by a complementary error function, multiplied with the Newtonian form. The resulting shortrange suppression compared to the Newtonian force is tabulated, because looking up from this table is faster than recomputing the corresponding factor, despite the memoryaccess panelty (which reduces cache performance) incurred by the table. Definition at line 1393 of file forcetree.c. References All, global_data_all_processes::ErrTolForceAcc, sph_particle_data::Hsml, particle_data::OldAcc, P, particle_data::Pos, SphP, and particle_data::Type. Referenced by gravity_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 2615 of file forcetree.c. Referenced by pmforce_periodic(), and pmpotential_periodic(). 

This function updates the toplevel tree after the multipole moments of the pseudoparticles have been updated. Definition at line 731 of file forcetree.c. References All, NODE::bitflags, NODE::d, DomainMoment, DomainNodeIndex, Extnodes, NODE::father, FLOAT, global_data_all_processes::ForceSoftening, NODE::mass, DomainNODE::mass, Nodes, NODE::s, DomainNODE::s, NODE::u, extNODE::vs, and DomainNODE::vs. 

This function updates the hmaxvalues in tree nodes that hold SPH particles. These values are needed to find all neighbors in the hydroforce computation. Since the Hsmlvalues are potentially changed in the SPHdenity computation, force_update_hmax() should be carried out just before the hydrodynamical SPH forces are computed, i.e. after density(). Definition at line 999 of file forcetree.c. References DomainHmax, DomainNodeIndex, and Extnodes. Referenced by compute_accelerations(). 

This function updates the sidelength of tree nodes in case the tree is not reconstructed, but only drifted. The grouping of particles to tree nodes is not changed in this case, but some tree nodes may need to be enlarged because particles moved out of their original bounds. Definition at line 870 of file forcetree.c. References DomainNodeIndex, DomainTreeNodeLen, and Nodes. Referenced by move_particles(). 

This routine updates the hmaxvalues of local tree nodes. Definition at line 1036 of file forcetree.c. References NODE::d, Extnodes, Father, NODE::father, extNODE::hmax, sph_particle_data::Hsml, Nodes, SphP, and NODE::u. 

This function recursively sets the hmaxvalues of the toplevel tree. Definition at line 1072 of file forcetree.c. References NODE::d, DomainHmax, DomainNodeIndex, Extnodes, NODE::father, extNODE::hmax, Nodes, and NODE::u. 

This function recursively enlarges nodes such that they always contain all their daughter nodes and daughter particles. Definition at line 908 of file forcetree.c. References NODE::center, NODE::d, Father, NODE::father, FLOAT, NODE::len, Nodes, P, particle_data::Pos, and NODE::u. 

This function recursively enlarges nodes of the toplevel tree such that they always contain all their daughter nodes. Definition at line 955 of file forcetree.c. References NODE::center, NODE::d, DomainNodeIndex, DomainTreeNodeLen, NODE::father, FLOAT, NODE::len, Nodes, and NODE::u. 

this routine determines the multipole moments for a given internal node and all its subnodes using a recursive computation. The result is stored in the Nodes[] structure in the sequence of this treewalk. Note that the bitflagsvariable for each node is used to store in the lowest bits some special information: Bit 0 flags whether the node belongs to the toplevel tree corresponding to the domain decomposition, while Bit 1 signals whether the toplevel node is dependent on local mass. If UNEQUALSOFTENINGS is set, bits 24 give the particle type with the maximum softening among the particles in the node, and bit 5 flags whether the node contains any particles with lower softening than that. Definition at line 425 of file forcetree.c. References All, NODE::bitflags, NODE::center, NODE::d, Extnodes, NODE::father, FLOAT, global_data_all_processes::ForceSoftening, extNODE::hmax, sph_particle_data::Hsml, NODE::mass, particle_data::Mass, global_data_all_processes::MaxPart, Nextnode, NODE::nextnode, Nodes, P, particle_data::Pos, NODE::s, NODE::sibling, SphP, NODE::suns, particle_data::Type, NODE::u, particle_data::Vel, and extNODE::vs. Referenced by force_treebuild_single(). 

This function updates the multipole moments of the pseudoparticles that represent the mass distribution on different CPUs. For that purpose, it first exchanges the necessary data, and then updates the toplevel tree accordingly. The detailed implementation of these two tasks is done in separate functions. Definition at line 674 of file forcetree.c. Referenced by force_treebuild(), and move_particles(). 