#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "allvars.h"
#include "proto.h"
#include "domain.h"
#include "forcetree.h"
#include "mymalloc.h"
#include "endrun.h"
/*! \file forcetree.c
* \brief gravitational tree and code for Ewald correction
*
* 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 domain 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. Tree nodes can
* be dynamically updated in drift/kick operations to avoid having to
* reconstruct the tree every timestep.
*/
struct NODE *Nodes_base, /*!< points to the actual memory allocted for the nodes */
*Nodes; /*!< this is a pointer used to access the nodes which is shifted such that Nodes[All.MaxPart]
gives the first allocated node */
struct extNODE *Extnodes, *Extnodes_base;
int MaxNodes; /*!< maximum allowed number of internal nodes */
int Numnodestree; /*!< number of (internal) nodes in each tree */
int *Nextnode; /*!< gives next node in tree walk (nodes array) */
int *Father; /*!< gives parent node in tree (Prenodes array) */
/*! auxiliary variable used to set-up non-recursive walk */
static int last;
static int tree_allocated_flag = 0;
static int
force_tree_build(int npart, struct unbind_data *mp);
static int
force_tree_build_single(int npart, struct unbind_data *mp);
static void
force_treeallocate(int maxnodes, int maxpart);
static void
force_update_hmax_of_node(int no, int mode);
static void
force_flag_localnodes(void);
static void
force_treeupdate_pseudos(int);
static void
force_update_node_recursive(int no, int sib, int father);
static void
force_create_empty_nodes(int no, int topnode, int bits, int x, int y, int z, int *nodecount, int *nextfree);
static void
force_exchange_pseudodata(void);
static void
force_insert_pseudo_particles(void);
int
force_tree_allocated()
{
return tree_allocated_flag;
}
void
force_tree_rebuild()
{
message(0, "Tree construction. (presently allocated=%g MB)\n", AllocatedBytes / (1024.0 * 1024.0));
if(force_tree_allocated()) {
force_tree_free();
}
/* construct tree if needed */
/* the tree is used in grav dens, hydro, bh and sfr */
force_treeallocate((int) (All.TreeAllocFactor * All.MaxPart) + NTopnodes, All.MaxPart);
walltime_measure("/Misc");
force_tree_build(NumPart, NULL);
walltime_measure("/Tree/Build");
message(0, "Tree construction done.\n");
}
/*! This function is a driver routine for constructing the gravitational
* oct-tree, which is done by calling a small number of other functions.
*/
int force_tree_build(int npart, struct unbind_data *mp)
{
int flag;
do
{
Numnodestree = force_tree_build_single(npart, mp);
MPI_Allreduce(&Numnodestree, &flag, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
if(flag == -1)
{
force_tree_free();
message(0, "Increasing TreeAllocFactor=%g", All.TreeAllocFactor);
All.TreeAllocFactor *= 1.15;
message(0, "new value=%g\n", All.TreeAllocFactor);
force_treeallocate((int) (All.TreeAllocFactor * All.MaxPart) + NTopnodes, All.MaxPart);
}
}
while(flag == -1);
force_flag_localnodes();
force_exchange_pseudodata();
force_treeupdate_pseudos(All.MaxPart);
TimeOfLastTreeConstruction = All.Time;
return Numnodestree;
}
/*! 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
* All.MaxPart.... All.MaxPart+nodes-1 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 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. 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.
*/
int force_tree_build_single(int npart, struct unbind_data *mp)
{
int i, j, k, subnode = 0, shift, parent, numnodes, rep;
int nfree, th, nn, no;
struct NODE *nfreep;
MyFloat lenhalf;
peanokey key, morton, th_key, *morton_list;
/* create an empty root node */
nfree = All.MaxPart; /* index of first free node */
nfreep = &Nodes[nfree]; /* select first node */
nfreep->len = DomainLen;
for(j = 0; j < 3; j++)
nfreep->center[j] = DomainCenter[j];
for(j = 0; j < 8; j++)
nfreep->u.suns[j] = -1;
numnodes = 1;
nfreep++;
nfree++;
/* create a set of empty nodes corresponding to the top-level domain
* 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 at the right place
*/
force_create_empty_nodes(All.MaxPart, 0, 1, 0, 0, 0, &numnodes, &nfree);
/* if a high-resolution region in a global tree is used, we need to generate
* an additional set empty nodes to make sure that we have a complete
* top-level tree for the high-resolution inset
*/
nfreep = &Nodes[nfree];
parent = -1; /* note: will not be used below before it is changed */
morton_list = (peanokey *) mymalloc("morton_list", NumPart * sizeof(peanokey));
/* now we insert all particles */
for(k = 0; k < npart; k++)
{
if(mp)
i = mp[k].index;
else
i = k;
/*We sometimes want to disable the tree for hot particles*/
if(P[i].Type == All.NoTreeType)
continue;
rep = 0;
key = peano_and_morton_key((int) ((P[i].Pos[0] - DomainCorner[0]) * DomainFac),
(int) ((P[i].Pos[1] - DomainCorner[1]) * DomainFac),
(int) ((P[i].Pos[2] - DomainCorner[2]) * DomainFac), BITS_PER_DIMENSION,
&morton);
morton_list[i] = morton;
shift = 3 * (BITS_PER_DIMENSION - 1);
no = 0;
while(TopNodes[no].Daughter >= 0)
{
no = TopNodes[no].Daughter + (key - TopNodes[no].StartKey) / (TopNodes[no].Size / 8);
shift -= 3;
}
no = TopNodes[no].Leaf;
th = DomainNodeIndex[no];
while(1)
{
if(th >= All.MaxPart) /* we are dealing with an internal node */
{
if(shift >= 0)
{
subnode = ((morton >> shift) & 7);
}
else
{
subnode = 0;
if(P[i].Pos[0] > Nodes[th].center[0])
subnode += 1;
if(P[i].Pos[1] > Nodes[th].center[1])
subnode += 2;
if(P[i].Pos[2] > Nodes[th].center[2])
subnode += 4;
}
#ifndef NOTREERND
if(Nodes[th].len < 1.0e-3 * All.ForceSoftening[P[i].Type])
{
/* seems like we're dealing with particles at identical (or extremely close)
* locations. Randomize subnode index to allow tree construction. Note: Multipole moments
* of tree are still correct, but this will only happen well below gravitational softening
* length-scale anyway.
*/
subnode = (int) (8.0 * get_random_number((P[i].ID + rep) % (RNDTABLE + (rep & 3))));
if(subnode >= 8)
subnode = 7;
rep++;
}
#endif
nn = Nodes[th].u.suns[subnode];
shift -= 3;
if(nn >= 0) /* ok, something is in the daughter slot already, need to continue */
{
parent = th;
th = nn;
}
else
{
/* here we have found an empty slot where we can attach
* the new particle as a leaf.
*/
Nodes[th].u.suns[subnode] = i;
break; /* done for this particle */
}
}
else
{
/* We try to insert into a leaf with a single particle. Need
* to generate a new internal node at this point.
*/
Nodes[parent].u.suns[subnode] = nfree;
nfreep->len = 0.5 * Nodes[parent].len;
lenhalf = 0.25 * Nodes[parent].len;
if(subnode & 1)
nfreep->center[0] = Nodes[parent].center[0] + lenhalf;
else
nfreep->center[0] = Nodes[parent].center[0] - lenhalf;
if(subnode & 2)
nfreep->center[1] = Nodes[parent].center[1] + lenhalf;
else
nfreep->center[1] = Nodes[parent].center[1] - lenhalf;
if(subnode & 4)
nfreep->center[2] = Nodes[parent].center[2] + lenhalf;
else
nfreep->center[2] = Nodes[parent].center[2] - lenhalf;
nfreep->u.suns[0] = -1;
nfreep->u.suns[1] = -1;
nfreep->u.suns[2] = -1;
nfreep->u.suns[3] = -1;
nfreep->u.suns[4] = -1;
nfreep->u.suns[5] = -1;
nfreep->u.suns[6] = -1;
nfreep->u.suns[7] = -1;
if(shift >= 0)
{
th_key = morton_list[th];
subnode = ((th_key >> shift) & 7);
}
else
{
subnode = 0;
if(P[th].Pos[0] > nfreep->center[0])
subnode += 1;
if(P[th].Pos[1] > nfreep->center[1])
subnode += 2;
if(P[th].Pos[2] > nfreep->center[2])
subnode += 4;
}
#ifndef NOTREERND
if(nfreep->len < 1.0e-3 * All.ForceSoftening[P[th].Type])
{
/* seems like we're dealing with particles at identical (or extremely close)
* locations. Randomize subnode index to allow tree construction. Note: Multipole moments
* of tree are still correct, but this will only happen well below gravitational softening
* length-scale anyway.
*/
subnode = (int) (8.0 * get_random_number((P[th].ID + rep) % (RNDTABLE + (rep & 3))));
if(subnode >= 8)
subnode = 7;
rep++;
}
#endif
nfreep->u.suns[subnode] = th;
th = nfree; /* resume trying to insert the new particle at
* the newly created internal node
*/
numnodes++;
nfree++;
nfreep++;
if((numnodes) >= MaxNodes)
{
message(1, "maximum number %d of tree-nodes reached for particle %d.\n",
MaxNodes, i);
if(All.TreeAllocFactor > 5.0)
{
message(1, "looks like a serious problem for particle %d, stopping with particle dump.\n", i);
savepositions(999999, 0);
endrun(1, "serious problem occured, snapshot saved.");
}
else
{
myfree(morton_list);
return -1;
}
}
}
}
}
myfree(morton_list);
/* insert the pseudo particles that represent the mass distribution of other domains */
force_insert_pseudo_particles();
/* now compute the multipole moments recursively */
last = -1;
force_update_node_recursive(All.MaxPart, -1, -1);
if(last >= All.MaxPart)
{
if(last >= All.MaxPart + MaxNodes) /* a pseudo-particle */
Nextnode[last - MaxNodes] = -1;
else
Nodes[last].u.d.nextnode = -1;
}
else
Nextnode[last] = -1;
return numnodes;
}
/*! This function recursively creates a set of empty tree nodes which
* corresponds to the top-level tree for the domain 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.
*/
void force_create_empty_nodes(int no, int topnode, int bits, int x, int y, int z, int *nodecount,
int *nextfree)
{
int i, j, k, n, sub, count;
MyFloat lenhalf;
if(TopNodes[topnode].Daughter >= 0)
{
for(i = 0; i < 2; i++)
for(j = 0; j < 2; j++)
for(k = 0; k < 2; k++)
{
sub = 7 & peano_hilbert_key((x << 1) + i, (y << 1) + j, (z << 1) + k, bits);
count = i + 2 * j + 4 * k;
Nodes[no].u.suns[count] = *nextfree;
lenhalf = 0.25 * Nodes[no].len;
Nodes[*nextfree].len = 0.5 * Nodes[no].len;
Nodes[*nextfree].center[0] = Nodes[no].center[0] + (2 * i - 1) * lenhalf;
Nodes[*nextfree].center[1] = Nodes[no].center[1] + (2 * j - 1) * lenhalf;
Nodes[*nextfree].center[2] = Nodes[no].center[2] + (2 * k - 1) * lenhalf;
for(n = 0; n < 8; n++)
Nodes[*nextfree].u.suns[n] = -1;
if(TopNodes[TopNodes[topnode].Daughter + sub].Daughter == -1)
DomainNodeIndex[TopNodes[TopNodes[topnode].Daughter + sub].Leaf] = *nextfree;
*nextfree = *nextfree + 1;
*nodecount = *nodecount + 1;
if((*nodecount) >= MaxNodes || (*nodecount) >= MaxTopNodes)
{
endrun(11, "maximum number MaxNodes=%d of tree-nodes reached."
"MaxTopNodes=%d NTopnodes=%d NTopleaves=%d nodecount=%d\n",
MaxNodes, MaxTopNodes, NTopnodes, NTopleaves, *nodecount);
}
force_create_empty_nodes(*nextfree - 1, TopNodes[topnode].Daughter + sub,
bits + 1, 2 * x + i, 2 * y + j, 2 * z + k, nodecount, nextfree);
}
}
}
/*! 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 domain-cell they correspond to. These quantities will be
* updated later on.
*/
void force_insert_pseudo_particles(void)
{
int i, index;
for(i = 0; i < NTopleaves; i++)
{
index = DomainNodeIndex[i];
if(DomainTask[i] != ThisTask)
Nodes[index].u.suns[0] = All.MaxPart + MaxNodes + i;
}
}
/*! 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 tree-walk.
*
* Note that the bitflags-variable for each node is used to store in the
* lowest bits some special information: Bit 0 flags whether the node
* belongs to the top-level tree corresponding to the domain
* decomposition, while Bit 1 signals whether the top-level node is
* dependent on local mass.
*
* bits 2-4 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.
*/
void force_update_node_recursive(int no, int sib, int father)
{
int j, jj, k, p, pp, nextsib, suns[8], count_particles, multiple_flag;
MyFloat hmax, vmax, v, divVmax;
MyFloat s[3], vs[3], mass;
struct particle_data *pa;
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
int maxsofttype, current_maxsofttype, diffsoftflag;
#else
MyFloat maxsoft;
#endif
if(no >= All.MaxPart && no < All.MaxPart + MaxNodes) /* internal node */
{
for(j = 0; j < 8; j++)
suns[j] = Nodes[no].u.suns[j]; /* this "backup" is necessary because the nextnode entry will overwrite one element (union!) */
if(last >= 0)
{
if(last >= All.MaxPart)
{
if(last >= All.MaxPart + MaxNodes) /* a pseudo-particle */
Nextnode[last - MaxNodes] = no;
else
Nodes[last].u.d.nextnode = no;
}
else
Nextnode[last] = no;
}
last = no;
mass = 0;
s[0] = 0;
s[1] = 0;
s[2] = 0;
vs[0] = 0;
vs[1] = 0;
vs[2] = 0;
hmax = 0;
vmax = 0;
divVmax = 0;
count_particles = 0;
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
maxsofttype = 7;
diffsoftflag = 0;
#else
maxsoft = 0;
#endif
for(j = 0; j < 8; j++)
{
if((p = suns[j]) >= 0)
{
/* check if we have a sibling on the same level */
for(jj = j + 1; jj < 8; jj++)
if((pp = suns[jj]) >= 0)
break;
if(jj < 8) /* yes, we do */
nextsib = pp;
else
nextsib = sib;
force_update_node_recursive(p, nextsib, no);
if(p >= All.MaxPart) /* an internal node or pseudo particle */
{
if(p >= All.MaxPart + MaxNodes) /* a pseudo particle */
{
/* nothing to be done here because the mass of the
* pseudo-particle is still zero. This will be changed
* later.
*/
}
else
{
mass += (Nodes[p].u.d.mass);
s[0] += (Nodes[p].u.d.mass * Nodes[p].u.d.s[0]);
s[1] += (Nodes[p].u.d.mass * Nodes[p].u.d.s[1]);
s[2] += (Nodes[p].u.d.mass * Nodes[p].u.d.s[2]);
vs[0] += (Nodes[p].u.d.mass * Extnodes[p].vs[0]);
vs[1] += (Nodes[p].u.d.mass * Extnodes[p].vs[1]);
vs[2] += (Nodes[p].u.d.mass * Extnodes[p].vs[2]);
if(Nodes[p].u.d.mass > 0)
{
if(Nodes[p].u.d.bitflags & (1 << BITFLAG_MULTIPLEPARTICLES))
count_particles += 2;
else
count_particles++;
}
if(Extnodes[p].hmax > hmax)
hmax = Extnodes[p].hmax;
if(Extnodes[p].vmax > vmax)
vmax = Extnodes[p].vmax;
if(Extnodes[p].divVmax > divVmax)
divVmax = Extnodes[p].divVmax;
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
diffsoftflag |= maskout_different_softening_flag(Nodes[p].u.d.bitflags);
if(maxsofttype == 7)
maxsofttype = extract_max_softening_type(Nodes[p].u.d.bitflags);
else
{
current_maxsofttype = extract_max_softening_type(Nodes[p].u.d.bitflags);
if(current_maxsofttype != 7)
{
if(All.ForceSoftening[current_maxsofttype] > All.ForceSoftening[maxsofttype])
{
maxsofttype = current_maxsofttype;
diffsoftflag = (1 << BITFLAG_MIXED_SOFTENINGS_IN_NODE);
}
else
{
if(All.ForceSoftening[current_maxsofttype] <
All.ForceSoftening[maxsofttype])
diffsoftflag = (1 << BITFLAG_MIXED_SOFTENINGS_IN_NODE);
}
}
}
#else
if(Nodes[p].maxsoft > maxsoft)
maxsoft = Nodes[p].maxsoft;
#endif
}
}
else /* a particle */
{
count_particles++;
pa = &P[p];
mass += (pa->Mass);
s[0] += (pa->Mass * pa->Pos[0]);
s[1] += (pa->Mass * pa->Pos[1]);
s[2] += (pa->Mass * pa->Pos[2]);
vs[0] += (pa->Mass * pa->Vel[0]);
vs[1] += (pa->Mass * pa->Vel[1]);
vs[2] += (pa->Mass * pa->Vel[2]);
if(pa->Type == 0)
{
if(P[p].Hsml > hmax)
hmax = P[p].Hsml;
if(SPHP(p).DivVel > divVmax)
divVmax = SPHP(p).DivVel;
}
for(k = 0; k < 3; k++)
if((v = fabs(pa->Vel[k])) > vmax)
vmax = v;
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
if(maxsofttype == 7)
maxsofttype = pa->Type;
else
{
if(All.ForceSoftening[pa->Type] > All.ForceSoftening[maxsofttype])
{
maxsofttype = pa->Type;
diffsoftflag = (1 << BITFLAG_MIXED_SOFTENINGS_IN_NODE);
}
else
{
if(All.ForceSoftening[pa->Type] < All.ForceSoftening[maxsofttype])
diffsoftflag = (1 << BITFLAG_MIXED_SOFTENINGS_IN_NODE);
}
}
#else
if(pa->Type == 0)
{
if(P[p].Hsml > maxsoft)
maxsoft = P[p].Hsml;
}
else
{
if(All.ForceSoftening[pa->Type] > maxsoft)
maxsoft = All.ForceSoftening[pa->Type];
}
#endif
}
}
}
if(mass)
{
s[0] /= mass;
s[1] /= mass;
s[2] /= mass;
vs[0] /= mass;
vs[1] /= mass;
vs[2] /= mass;
}
else
{
s[0] = Nodes[no].center[0];
s[1] = Nodes[no].center[1];
s[2] = Nodes[no].center[2];
vs[0] = 0;
vs[1] = 0;
vs[2] = 0;
}
Nodes[no].Ti_current = All.Ti_Current;
Nodes[no].u.d.mass = mass;
Nodes[no].u.d.s[0] = s[0];
Nodes[no].u.d.s[1] = s[1];
Nodes[no].u.d.s[2] = s[2];
Extnodes[no].Ti_lastkicked = All.Ti_Current;
Extnodes[no].Flag = GlobFlag;
Extnodes[no].vs[0] = vs[0];
Extnodes[no].vs[1] = vs[1];
Extnodes[no].vs[2] = vs[2];
Extnodes[no].hmax = hmax;
Extnodes[no].vmax = vmax;
Extnodes[no].divVmax = divVmax;
Extnodes[no].dp[0] = 0;
Extnodes[no].dp[1] = 0;
Extnodes[no].dp[2] = 0;
if(count_particles > 1) /* this flags that the node represents more than one particle */
multiple_flag = (1 << BITFLAG_MULTIPLEPARTICLES);
else
multiple_flag = 0;
Nodes[no].u.d.bitflags = multiple_flag;
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
Nodes[no].u.d.bitflags |= diffsoftflag + (maxsofttype << BITFLAG_MAX_SOFTENING_TYPE);
#else
Nodes[no].maxsoft = maxsoft;
#endif
Nodes[no].u.d.sibling = sib;
Nodes[no].u.d.father = father;
}
else /* single particle or pseudo particle */
{
if(last >= 0)
{
if(last >= All.MaxPart)
{
if(last >= All.MaxPart + MaxNodes) /* a pseudo-particle */
Nextnode[last - MaxNodes] = no;
else
Nodes[last].u.d.nextnode = no;
}
else
Nextnode[last] = no;
}
last = no;
if(no < All.MaxPart) /* only set it for single particles */
Father[no] = father;
}
}
/*! This function communicates the values of the multipole moments of the
* top-level tree-nodes of the domain grid. This data can then be used to
* update the pseudo-particles on each CPU accordingly.
*/
void force_exchange_pseudodata(void)
{
int i, no, m, ta, recvTask;
int *recvcounts, *recvoffset;
struct DomainNODE
{
MyFloat s[3];
MyFloat vs[3];
MyFloat mass;
MyFloat hmax;
MyFloat vmax;
MyFloat divVmax;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
MyFloat maxsoft;
#endif
unsigned int bitflags;
}
*DomainMoment;
DomainMoment = (struct DomainNODE *) mymalloc("DomainMoment", NTopleaves * sizeof(struct DomainNODE));
memset(&DomainMoment[0], 0, sizeof(DomainMoment[0]) * NTopleaves);
for(m = 0; m < All.DomainOverDecompositionFactor; m++)
for(i = DomainStartList[ThisTask * All.DomainOverDecompositionFactor + m];
i <= DomainEndList[ThisTask * All.DomainOverDecompositionFactor + m]; i++)
{
no = DomainNodeIndex[i];
/* read out the multipole moments from the local base cells */
DomainMoment[i].s[0] = Nodes[no].u.d.s[0];
DomainMoment[i].s[1] = Nodes[no].u.d.s[1];
DomainMoment[i].s[2] = Nodes[no].u.d.s[2];
DomainMoment[i].vs[0] = Extnodes[no].vs[0];
DomainMoment[i].vs[1] = Extnodes[no].vs[1];
DomainMoment[i].vs[2] = Extnodes[no].vs[2];
DomainMoment[i].mass = Nodes[no].u.d.mass;
DomainMoment[i].hmax = Extnodes[no].hmax;
DomainMoment[i].vmax = Extnodes[no].vmax;
DomainMoment[i].divVmax = Extnodes[no].divVmax;
DomainMoment[i].bitflags = Nodes[no].u.d.bitflags;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
DomainMoment[i].maxsoft = Nodes[no].maxsoft;
#endif
}
/* share the pseudo-particle data accross CPUs */
recvcounts = (int *) mymalloc("recvcounts", sizeof(int) * NTask);
recvoffset = (int *) mymalloc("recvoffset", sizeof(int) * NTask);
for(m = 0; m < All.DomainOverDecompositionFactor; m++)
{
for(recvTask = 0; recvTask < NTask; recvTask++)
{
recvcounts[recvTask] =
(DomainEndList[recvTask * All.DomainOverDecompositionFactor + m]
- DomainStartList[recvTask * All.DomainOverDecompositionFactor + m] +
1)
* sizeof(struct DomainNODE);
recvoffset[recvTask] = DomainStartList[recvTask * All.DomainOverDecompositionFactor + m] * sizeof(struct DomainNODE);
}
MPI_Allgatherv(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL,
&DomainMoment[0], recvcounts, recvoffset,
MPI_BYTE, MPI_COMM_WORLD);
}
myfree(recvoffset);
myfree(recvcounts);
for(ta = 0; ta < NTask; ta++)
if(ta != ThisTask)
for(m = 0; m < All.DomainOverDecompositionFactor; m++)
for(i = DomainStartList[ta * All.DomainOverDecompositionFactor + m]; i <= DomainEndList[ta * All.DomainOverDecompositionFactor + m]; i++)
{
no = DomainNodeIndex[i];
Nodes[no].u.d.s[0] = DomainMoment[i].s[0];
Nodes[no].u.d.s[1] = DomainMoment[i].s[1];
Nodes[no].u.d.s[2] = DomainMoment[i].s[2];
Extnodes[no].vs[0] = DomainMoment[i].vs[0];
Extnodes[no].vs[1] = DomainMoment[i].vs[1];
Extnodes[no].vs[2] = DomainMoment[i].vs[2];
Nodes[no].u.d.mass = DomainMoment[i].mass;
Extnodes[no].hmax = DomainMoment[i].hmax;
Extnodes[no].vmax = DomainMoment[i].vmax;
Extnodes[no].divVmax = DomainMoment[i].divVmax;
Nodes[no].u.d.bitflags =
(Nodes[no].u.d.bitflags & (~BITFLAG_MASK)) | (DomainMoment[i].bitflags & BITFLAG_MASK);
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
Nodes[no].maxsoft = DomainMoment[i].maxsoft;
#endif
}
myfree(DomainMoment);
}
/*! This function updates the top-level tree after the multipole moments of
* the pseudo-particles have been updated.
*/
void force_treeupdate_pseudos(int no)
{
int j, p, count_particles, multiple_flag;
MyFloat hmax, vmax, divVmax;
MyFloat s[3], vs[3], mass;
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
int maxsofttype, diffsoftflag, current_maxsofttype;
#else
MyFloat maxsoft;
#endif
mass = 0;
s[0] = 0;
s[1] = 0;
s[2] = 0;
vs[0] = 0;
vs[1] = 0;
vs[2] = 0;
hmax = 0;
vmax = 0;
divVmax = 0;
count_particles = 0;
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
maxsofttype = 7;
diffsoftflag = 0;
#else
maxsoft = 0;
#endif
p = Nodes[no].u.d.nextnode;
for(j = 0; j < 8; j++) /* since we are dealing with top-level nodes, we now that there are 8 consecutive daughter nodes */
{
if(p >= All.MaxPart && p < All.MaxPart + MaxNodes) /* internal node */
{
if(Nodes[p].u.d.bitflags & (1 << BITFLAG_INTERNAL_TOPLEVEL))
force_treeupdate_pseudos(p);
mass += (Nodes[p].u.d.mass);
s[0] += (Nodes[p].u.d.mass * Nodes[p].u.d.s[0]);
s[1] += (Nodes[p].u.d.mass * Nodes[p].u.d.s[1]);
s[2] += (Nodes[p].u.d.mass * Nodes[p].u.d.s[2]);
vs[0] += (Nodes[p].u.d.mass * Extnodes[p].vs[0]);
vs[1] += (Nodes[p].u.d.mass * Extnodes[p].vs[1]);
vs[2] += (Nodes[p].u.d.mass * Extnodes[p].vs[2]);
if(Extnodes[p].hmax > hmax)
hmax = Extnodes[p].hmax;
if(Extnodes[p].vmax > vmax)
vmax = Extnodes[p].vmax;
if(Extnodes[p].divVmax > divVmax)
divVmax = Extnodes[p].divVmax;
if(Nodes[p].u.d.mass > 0)
{
if(Nodes[p].u.d.bitflags & (1 << BITFLAG_MULTIPLEPARTICLES))
count_particles += 2;
else
count_particles++;
}
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
diffsoftflag |= maskout_different_softening_flag(Nodes[p].u.d.bitflags);
if(maxsofttype == 7)
maxsofttype = extract_max_softening_type(Nodes[p].u.d.bitflags);
else
{
current_maxsofttype = extract_max_softening_type(Nodes[p].u.d.bitflags);
if(current_maxsofttype != 7)
{
if(All.ForceSoftening[current_maxsofttype] > All.ForceSoftening[maxsofttype])
{
maxsofttype = current_maxsofttype;
diffsoftflag = (1 << BITFLAG_MIXED_SOFTENINGS_IN_NODE);
}
else
{
if(All.ForceSoftening[current_maxsofttype] < All.ForceSoftening[maxsofttype])
diffsoftflag = (1 << BITFLAG_MIXED_SOFTENINGS_IN_NODE);
}
}
}
#else
if(Nodes[p].maxsoft > maxsoft)
maxsoft = Nodes[p].maxsoft;
#endif
}
else
endrun(6767, "may not happen"); /* may not happen */
p = Nodes[p].u.d.sibling;
}
if(mass)
{
s[0] /= mass;
s[1] /= mass;
s[2] /= mass;
vs[0] /= mass;
vs[1] /= mass;
vs[2] /= mass;
}
else
{
s[0] = Nodes[no].center[0];
s[1] = Nodes[no].center[1];
s[2] = Nodes[no].center[2];
vs[0] = 0;
vs[1] = 0;
vs[2] = 0;
}
Nodes[no].u.d.s[0] = s[0];
Nodes[no].u.d.s[1] = s[1];
Nodes[no].u.d.s[2] = s[2];
Extnodes[no].vs[0] = vs[0];
Extnodes[no].vs[1] = vs[1];
Extnodes[no].vs[2] = vs[2];
Nodes[no].u.d.mass = mass;
Extnodes[no].hmax = hmax;
Extnodes[no].vmax = vmax;
Extnodes[no].divVmax = divVmax;
Extnodes[no].Flag = GlobFlag;
if(count_particles > 1)
multiple_flag = (1 << BITFLAG_MULTIPLEPARTICLES);
else
multiple_flag = 0;
Nodes[no].u.d.bitflags &= (~BITFLAG_MASK); /* this clears the bits */
Nodes[no].u.d.bitflags |= multiple_flag;
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
Nodes[no].u.d.bitflags |= diffsoftflag + (maxsofttype << BITFLAG_MAX_SOFTENING_TYPE);
#else
Nodes[no].maxsoft = maxsoft;
#endif
}
/*! This function flags nodes in the top-level tree that are dependent on
* local particle data.
*/
static void
force_flag_localnodes(void)
{
int no, i, m;
/* mark all top-level nodes */
for(i = 0; i < NTopleaves; i++)
{
no = DomainNodeIndex[i];
while(no >= 0)
{
if(Nodes[no].u.d.bitflags & (1 << BITFLAG_TOPLEVEL))
break;
Nodes[no].u.d.bitflags |= (1 << BITFLAG_TOPLEVEL);
no = Nodes[no].u.d.father;
}
/* mark also internal top level nodes */
no = DomainNodeIndex[i];
no = Nodes[no].u.d.father;
while(no >= 0)
{
if(Nodes[no].u.d.bitflags & (1 << BITFLAG_INTERNAL_TOPLEVEL))
break;
Nodes[no].u.d.bitflags |= (1 << BITFLAG_INTERNAL_TOPLEVEL);
no = Nodes[no].u.d.father;
}
}
/* mark top-level nodes that contain local particles */
for(m = 0; m < All.DomainOverDecompositionFactor; m++)
for(i = DomainStartList[ThisTask * All.DomainOverDecompositionFactor + m];
i <= DomainEndList[ThisTask * All.DomainOverDecompositionFactor + m]; i++)
{
no = DomainNodeIndex[i];
if(DomainTask[i] != ThisTask)
endrun(131231231, "DomainTask struct is corrupted");
while(no >= 0)
{
if(Nodes[no].u.d.bitflags & (1 << BITFLAG_DEPENDS_ON_LOCAL_MASS))
break;
Nodes[no].u.d.bitflags |= (1 << BITFLAG_DEPENDS_ON_LOCAL_MASS);
no = Nodes[no].u.d.father;
}
}
}
static void real_force_drift_node(int no, int time1);
void force_drift_node(int no, int time1) {
force_drift_node_full(no, time1, 1);
}
int force_drift_node_full(int no, int time1, int blocking) {
if(time1 == Nodes[no].Ti_current)
return 0;
#pragma omp atomic
TotalNodeDrifts ++;
#ifdef OPENMP_USE_SPINLOCK
int lockstate;
if (blocking) {
lockstate = pthread_spin_lock(&Nodes[no].SpinLock);
} else {
lockstate = pthread_spin_trylock(&Nodes[no].SpinLock);
}
if(0 == lockstate) {
if(time1 != Nodes[no].Ti_current) {
real_force_drift_node(no, time1);
#pragma omp flush
} else {
#pragma omp atomic
BlockedNodeDrifts ++;
}
pthread_spin_unlock(&Nodes[no].SpinLock);
return 0;
} else {
if(blocking) {
endrun(99999, "shall not happend");
}
return -1;
}
#else
/* do not use spinlock */
#pragma omp critical (_driftnode_)
{
if(time1 != Nodes[no].Ti_current) {
real_force_drift_node(no, time1);
} else {
BlockedNodeDrifts ++;
}
}
return 0;
#endif
}
static void real_force_drift_node(int no, int time1)
{
int j;
double dt_drift, dt_drift_hmax, fac;
if(time1 == Nodes[no].Ti_current) return;
if(Nodes[no].u.d.bitflags & (1 << BITFLAG_NODEHASBEENKICKED))
{
if(Extnodes[no].Ti_lastkicked != Nodes[no].Ti_current)
{
endrun(1, "inconsistency in drift node: Extnodes[no].Ti_lastkicked=%d Nodes[no].Ti_current=%d\n",
Extnodes[no].Ti_lastkicked, Nodes[no].Ti_current);
}
if(Nodes[no].u.d.mass)
fac = 1 / Nodes[no].u.d.mass;
else
fac = 0;
for(j = 0; j < 3; j++)
{
Extnodes[no].vs[j] += fac * (Extnodes[no].dp[j]);
Extnodes[no].dp[j] = 0;
}
Nodes[no].u.d.bitflags &= (~(1 << BITFLAG_NODEHASBEENKICKED));
}
dt_drift_hmax = get_drift_factor(Nodes[no].Ti_current, time1);
dt_drift = dt_drift_hmax;
for(j = 0; j < 3; j++)
Nodes[no].u.d.s[j] += Extnodes[no].vs[j] * dt_drift;
Nodes[no].len += 2 * Extnodes[no].vmax * dt_drift;
// Extnodes[no].hmax *= exp(0.333333333333 * Extnodes[no].divVmax * dt_drift_hmax);
Nodes[no].Ti_current = time1;
}
void force_kick_node(int i, MyFloat * dv)
{
int j, no;
MyFloat dp[3], v, vmax;
/*We sometimes want to disable the tree for hot particles*/
if(P[i].Type == All.NoTreeType)
return;
for(j = 0; j < 3; j++)
{
dp[j] = P[i].Mass * dv[j];
}
for(j = 0, vmax = 0; j < 3; j++)
if((v = fabs(P[i].Vel[j])) > vmax)
vmax = v;
no = Father[i];
while(no >= 0)
{
real_force_drift_node(no, All.Ti_Current);
for(j = 0; j < 3; j++)
{
Extnodes[no].dp[j] += dp[j];
}
if(Extnodes[no].vmax < vmax)
Extnodes[no].vmax = vmax;
Nodes[no].u.d.bitflags |= (1 << BITFLAG_NODEHASBEENKICKED);
Extnodes[no].Ti_lastkicked = All.Ti_Current;
if(Nodes[no].u.d.bitflags & (1 << BITFLAG_TOPLEVEL)) /* top-level tree-node reached */
{
if(Extnodes[no].Flag != GlobFlag)
{
Extnodes[no].Flag = GlobFlag;
DomainList[DomainNumChanged++] = no;
}
break;
}
no = Nodes[no].u.d.father;
}
}
void force_finish_kick_nodes(void)
{
int i, j, no, ta, totDomainNumChanged;
int *domainList_all;
int *counts, *counts_dp, *offset_list, *offset_dp, *offset_vmax;
MyDouble *domainDp_loc, *domainDp_all;
MyFloat *domainVmax_loc, *domainVmax_all;
/* share the momentum-data of the pseudo-particles accross CPUs */
counts = (int *) mymalloc("counts", sizeof(int) * NTask);
counts_dp = (int *) mymalloc("counts_dp", sizeof(int) * NTask);
offset_list = (int *) mymalloc("offset_list", sizeof(int) * NTask);
offset_dp = (int *) mymalloc("offset_dp", sizeof(int) * NTask);
offset_vmax = (int *) mymalloc("offset_vmax", sizeof(int) * NTask);
domainDp_loc = (MyDouble *) mymalloc("domainDp_loc", DomainNumChanged * 3 * sizeof(MyDouble));
domainVmax_loc = (MyFloat *) mymalloc("domainVmax_loc", DomainNumChanged * sizeof(MyFloat));
for(i = 0; i < DomainNumChanged; i++)
{
for(j = 0; j < 3; j++)
{
domainDp_loc[i * 3 + j] = Extnodes[DomainList[i]].dp[j];
}
domainVmax_loc[i] = Extnodes[DomainList[i]].vmax;
}
MPI_Allgather(&DomainNumChanged, 1, MPI_INT, counts, 1, MPI_INT, MPI_COMM_WORLD);
for(ta = 0, totDomainNumChanged = 0, offset_list[0] = 0, offset_dp[0] = 0, offset_vmax[0] = 0; ta < NTask;
ta++)
{
totDomainNumChanged += counts[ta];
if(ta > 0)
{
offset_list[ta] = offset_list[ta - 1] + counts[ta - 1];
offset_dp[ta] = offset_dp[ta - 1] + counts[ta - 1] * 3 * sizeof(MyDouble);
offset_vmax[ta] = offset_vmax[ta - 1] + counts[ta - 1] * sizeof(MyFloat);
}
}
message(0, "I exchange kick momenta for %d top-level nodes out of %d\n", totDomainNumChanged, NTopleaves);
domainDp_all = (MyDouble *) mymalloc("domainDp_all", totDomainNumChanged * 3 * sizeof(MyDouble));
domainVmax_all = (MyFloat *) mymalloc("domainVmax_all", totDomainNumChanged * sizeof(MyFloat));
domainList_all = (int *) mymalloc("domainList_all", totDomainNumChanged * sizeof(int));
MPI_Allgatherv(DomainList, DomainNumChanged, MPI_INT,
domainList_all, counts, offset_list, MPI_INT, MPI_COMM_WORLD);
for(ta = 0; ta < NTask; ta++)
{
counts_dp[ta] = counts[ta] * 3 * sizeof(MyDouble);
counts[ta] *= sizeof(MyFloat);
}
MPI_Allgatherv(domainDp_loc, DomainNumChanged * 3 * sizeof(MyDouble), MPI_BYTE,
domainDp_all, counts_dp, offset_dp, MPI_BYTE, MPI_COMM_WORLD);
MPI_Allgatherv(domainVmax_loc, DomainNumChanged * sizeof(MyFloat), MPI_BYTE,
domainVmax_all, counts, offset_vmax, MPI_BYTE, MPI_COMM_WORLD);
/* construct momentum kicks in top-level tree */
for(i = 0; i < totDomainNumChanged; i++)
{
no = domainList_all[i];
if(Nodes[no].u.d.bitflags & (1 << BITFLAG_DEPENDS_ON_LOCAL_MASS)) /* to avoid that the local one is kicked twice */
no = Nodes[no].u.d.father;
while(no >= 0)
{
real_force_drift_node(no, All.Ti_Current);
for(j = 0; j < 3; j++)
{
Extnodes[no].dp[j] += domainDp_all[3 * i + j];
}
if(Extnodes[no].vmax < domainVmax_all[i])
Extnodes[no].vmax = domainVmax_all[i];
Nodes[no].u.d.bitflags |= (1 << BITFLAG_NODEHASBEENKICKED);
Extnodes[no].Ti_lastkicked = All.Ti_Current;
no = Nodes[no].u.d.father;
}
}
myfree(domainList_all);
myfree(domainVmax_all);
myfree(domainDp_all);
myfree(domainVmax_loc);
myfree(domainDp_loc);
myfree(offset_vmax);
myfree(offset_dp);
myfree(offset_list);
myfree(counts_dp);
myfree(counts);
}
/*! This function updates the hmax-values in tree nodes that hold SPH
* particles. These values are needed to find all neighbors in the
* hydro-force computation. Since the Hsml-values are potentially changed
* in the SPH-denity computation, force_update_hmax() should be carried
* out just before the hydrodynamical SPH forces are computed, i.e. after
* density().
*/
void force_update_hmax(void)
{
int i, no, ta, totDomainNumChanged;
int *domainList_all;
int *counts, *offset_list, *offset_hmax;
MyFloat *domainHmax_loc, *domainHmax_all;
walltime_measure("/Misc");
GlobFlag++;
DomainNumChanged = 0;
DomainList = (int *) mymalloc("DomainList", NTopleaves * sizeof(int));
for(i = FirstActiveParticle; i >= 0; i = NextActiveParticle[i])
if(P[i].Type == 0)
{
no = Father[i];
while(no >= 0)
{
real_force_drift_node(no, All.Ti_Current);
if(P[i].Hsml > Extnodes[no].hmax || SPHP(i).DivVel > Extnodes[no].divVmax)
{
if(P[i].Hsml > Extnodes[no].hmax)
Extnodes[no].hmax = P[i].Hsml;
if(SPHP(i).DivVel > Extnodes[no].divVmax)
Extnodes[no].divVmax = SPHP(i).DivVel;
if(Nodes[no].u.d.bitflags & (1 << BITFLAG_TOPLEVEL)) /* we reached a top-level node */
{
if(Extnodes[no].Flag != GlobFlag)
{
Extnodes[no].Flag = GlobFlag;
DomainList[DomainNumChanged++] = no;
}
break;
}
}
else
break;
no = Nodes[no].u.d.father;
}
}
/* share the hmax-data of the pseudo-particles accross CPUs */
counts = (int *) mymalloc("counts", sizeof(int) * NTask);
offset_list = (int *) mymalloc("offset_list", sizeof(int) * NTask);
offset_hmax = (int *) mymalloc("offset_hmax", sizeof(int) * NTask);
domainHmax_loc = (MyFloat *) mymalloc("domainHmax_loc", DomainNumChanged * 2 * sizeof(MyFloat));
for(i = 0; i < DomainNumChanged; i++)
{
domainHmax_loc[2 * i] = Extnodes[DomainList[i]].hmax;
domainHmax_loc[2 * i + 1] = Extnodes[DomainList[i]].divVmax;
}
MPI_Allgather(&DomainNumChanged, 1, MPI_INT, counts, 1, MPI_INT, MPI_COMM_WORLD);
for(ta = 0, totDomainNumChanged = 0, offset_list[0] = 0, offset_hmax[0] = 0; ta < NTask; ta++)
{
totDomainNumChanged += counts[ta];
if(ta > 0)
{
offset_list[ta] = offset_list[ta - 1] + counts[ta - 1];
offset_hmax[ta] = offset_hmax[ta - 1] + counts[ta - 1] * 2 * sizeof(MyFloat);
}
}
message(0, "Hmax exchange: %d topleaves out of %d\n", totDomainNumChanged, NTopleaves);
domainHmax_all = (MyFloat *) mymalloc("domainHmax_all", totDomainNumChanged * 2 * sizeof(MyFloat));
domainList_all = (int *) mymalloc("domainList_all", totDomainNumChanged * sizeof(int));
MPI_Allgatherv(DomainList, DomainNumChanged, MPI_INT,
domainList_all, counts, offset_list, MPI_INT, MPI_COMM_WORLD);
for(ta = 0; ta < NTask; ta++)
counts[ta] *= 2 * sizeof(MyFloat);
MPI_Allgatherv(domainHmax_loc, 2 * DomainNumChanged * sizeof(MyFloat), MPI_BYTE,
domainHmax_all, counts, offset_hmax, MPI_BYTE, MPI_COMM_WORLD);
for(i = 0; i < totDomainNumChanged; i++)
{
no = domainList_all[i];
if(Nodes[no].u.d.bitflags & (1 << BITFLAG_DEPENDS_ON_LOCAL_MASS)) /* to avoid that the hmax is updated twice */
no = Nodes[no].u.d.father;
while(no >= 0)
{
real_force_drift_node(no, All.Ti_Current);
if(domainHmax_all[2 * i] > Extnodes[no].hmax || domainHmax_all[2 * i + 1] > Extnodes[no].divVmax)
{
if(domainHmax_all[2 * i] > Extnodes[no].hmax)
Extnodes[no].hmax = domainHmax_all[2 * i];
if(domainHmax_all[2 * i + 1] > Extnodes[no].divVmax)
Extnodes[no].divVmax = domainHmax_all[2 * i + 1];
}
else
break;
no = Nodes[no].u.d.father;
}
}
myfree(domainList_all);
myfree(domainHmax_all);
myfree(domainHmax_loc);
myfree(offset_hmax);
myfree(offset_list);
myfree(counts);
myfree(DomainList);
walltime_measure("/Tree/HmaxUpdate");
}
/*! 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.
*/
void force_treeallocate(int maxnodes, int maxpart)
{
size_t bytes;
double allbytes = 0, allbytes_topleaves = 0;
tree_allocated_flag = 1;
DomainNodeIndex = (int *) mymalloc("DomainNodeIndex", bytes = NTopleaves * sizeof(int));
allbytes_topleaves += bytes;
MaxNodes = maxnodes;
if(!(Nodes_base = (struct NODE *) mymalloc("Nodes_base", bytes = (MaxNodes + 1) * sizeof(struct NODE))))
{
endrun(3, "failed to allocate memory for %d tree-nodes (%g MB).\n", MaxNodes, bytes / (1024.0 * 1024.0));
}
allbytes += bytes;
if(!
(Extnodes_base =
(struct extNODE *) mymalloc("Extnodes_base", bytes = (MaxNodes + 1) * sizeof(struct extNODE))))
{
endrun(3, "failed to allocate memory for %d tree-extnodes (%g MB).\n",
MaxNodes, bytes / (1024.0 * 1024.0));
}
#ifdef OPENMP_USE_SPINLOCK
{
int i;
for (i = 0; i < MaxNodes + 1; i ++) {
/* Maybe we can directly set these guys to one
*
* at least with the glibc spinlock implementation.
* */
pthread_spin_init(&Nodes_base[i].SpinLock, 0);
}
}
#endif
allbytes += bytes;
Nodes = Nodes_base - All.MaxPart;
Extnodes = Extnodes_base - All.MaxPart;
if(!(Nextnode = (int *) mymalloc("Nextnode", bytes = (maxpart + NTopnodes) * sizeof(int))))
{
endrun(1, "Failed to allocate %d spaces for 'Nextnode' array (%g MB)\n",
maxpart + NTopnodes, bytes / (1024.0 * 1024.0));
}
allbytes += bytes;
if(!(Father = (int *) mymalloc("Father", bytes = (maxpart) * sizeof(int))))
{
endrun(1, "Failed to allocate %d spaces for 'Father' array (%g MB)\n", maxpart, bytes / (1024.0 * 1024.0));
}
allbytes += bytes;
message(0, "Allocated %g MByte for BH-tree, and %g Mbyte for top-leaves. (presently allocted %g MB)\n",
allbytes / (1024.0 * 1024.0), allbytes_topleaves / (1024.0 * 1024.0),
AllocatedBytes / (1024.0 * 1024.0));
}
/*! This function frees the memory allocated for the tree, i.e. it frees
* the space allocated by the function force_treeallocate().
*/
void force_tree_free(void)
{
myfree(Father);
myfree(Nextnode);
myfree(Extnodes_base);
myfree(Nodes_base);
myfree(DomainNodeIndex);
tree_allocated_flag = 0;
}