Homec4science
Diffusion Gear (master) Browse

Gear/src/forcetree_sub.c
61b099c90e56master

Last Change
Raw File

forcetree_sub.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
#ifdef PY_INTERFACE
/*! \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.
*/
/*! auxialiary variable used to set-up non-recursive walk */
static int last;
/*! length of lock-up table for short-range force kernel in TreePM algorithm */
#define NTAB 1000
/*! variables for short-range lookup table */
static float tabfac, shortrange_table[NTAB], shortrange_table_potential[NTAB];
/*! toggles after first tree-memory allocation, has only influence on log-files */
static int first_flag = 0;
#ifdef PERIODIC
/*! Macro that maps a distance to the nearest periodic neighbour */
#define NEAREST(x) (((x)>boxhalf)?((x)-boxsize):(((x)<-boxhalf)?((x)+boxsize):(x)))
/*! Size of 3D lock-up table for Ewald correction force */
#define EN 64
/*! 3D lock-up table for Ewald correction to force and potential. Only one
* octant is stored, the rest constructed by using the symmetry
*/
static FLOAT fcorrx[EN + 1][EN + 1][EN + 1];
static FLOAT fcorry[EN + 1][EN + 1][EN + 1];
static FLOAT fcorrz[EN + 1][EN + 1][EN + 1];
static FLOAT potcorr[EN + 1][EN + 1][EN + 1];
static double fac_intp;
#endif
/*! 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
* cell-opening criterion, or the `relative' opening criterion is used.
*/
int force_treeevaluate_sub(int target, int mode, double *ewaldcountsum)
{
struct NODE *nop = 0;
int no, ninteractions, ptype;
double r2, dx, dy, dz, mass, r, fac, u, h, h_inv, h3_inv;
double acc_x, acc_y, acc_z, pos_x, pos_y, pos_z, aold;
#if defined(UNEQUALSOFTENINGS) && !defined(ADAPTIVE_GRAVSOFT_FORGAS)
int maxsofttype;
#endif
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
double soft = 0;
#endif
#ifdef PERIODIC
double boxsize, boxhalf;
boxsize = All.BoxSize;
boxhalf = 0.5 * All.BoxSize;
#endif
acc_x = 0;
acc_y = 0;
acc_z = 0;
ninteractions = 0;
if(mode == 0)
{
pos_x = Q[target].Pos[0];
pos_y = Q[target].Pos[1];
pos_z = Q[target].Pos[2];
ptype = Q[target].Type;
aold = All.ErrTolForceAcc * Q[target].OldAcc;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
soft = SphQ[target].Hsml;
#endif
}
else
{
pos_x = GravDataGet[target].u.Pos[0];
pos_y = GravDataGet[target].u.Pos[1];
pos_z = GravDataGet[target].u.Pos[2];
#ifdef UNEQUALSOFTENINGS
ptype = GravDataGet[target].Type;
#else
ptype = P[0].Type;
#endif
aold = All.ErrTolForceAcc * GravDataGet[target].w.OldAcc;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
soft = GravDataGet[target].Soft;
#endif
}
#ifndef UNEQUALSOFTENINGS
h = All.ForceSofteningQ;
h_inv = 1.0 / h;
h3_inv = h_inv * h_inv * h_inv;
#endif
no = All.MaxPart; /* root node */
while(no >= 0)
{
if(no < All.MaxPart) /* single particle */
{
/* the index of the node is the index of the particle */
/* observe the sign */
dx = P[no].Pos[0] - pos_x;
dy = P[no].Pos[1] - pos_y;
dz = P[no].Pos[2] - pos_z;
mass = P[no].Mass;
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
if(mode == 0)
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
}
no = Nextnode[no - MaxNodes];
continue;
}
nop = &Nodes[no];
dx = nop->u.d.s[0] - pos_x;
dy = nop->u.d.s[1] - pos_y;
dz = nop->u.d.s[2] - pos_z;
mass = nop->u.d.mass;
}
#ifdef PERIODIC
dx = NEAREST(dx);
dy = NEAREST(dy);
dz = NEAREST(dz);
#endif
r2 = dx * dx + dy * dy + dz * dz;
if(no < All.MaxPart)
{
#ifdef UNEQUALSOFTENINGS
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
h = soft;
else
h = All.ForceSofteningQ;
if(P[no].Type == 0)
{
if(h < SphP[no].Hsml)
h = SphP[no].Hsml;
}
else
{
if(h < All.ForceSoftening[P[no].Type])
h = All.ForceSoftening[P[no].Type];
}
#else
h = All.ForceSofteningQ;
if(h < All.ForceSoftening[P[no].Type])
h = All.ForceSoftening[P[no].Type];
#endif
#endif
no = Nextnode[no];
}
else /* we have an internal node. Need to check opening criterion */
{
if(mode == 1)
{
if((nop->u.d.bitflags & 3) == 1) /* if it's a top-level node
* which does not contain
* local particles we can
* continue to do a short-cut */
{
no = nop->u.d.sibling;
continue;
}
}
if(All.ErrTolTheta) /* check Barnes-Hut opening criterion */
{
if(nop->len * nop->len > r2 * All.ErrTolTheta * All.ErrTolTheta)
{
/* open cell */
no = nop->u.d.nextnode;
continue;
}
}
else /* check relative opening criterion */
{
if(mass * nop->len * nop->len > r2 * r2 * aold)
{
/* open cell */
no = nop->u.d.nextnode;
continue;
}
/* check in addition whether we lie inside the cell */
if(fabs(nop->center[0] - pos_x) < 0.60 * nop->len)
{
if(fabs(nop->center[1] - pos_y) < 0.60 * nop->len)
{
if(fabs(nop->center[2] - pos_z) < 0.60 * nop->len)
{
no = nop->u.d.nextnode;
continue;
}
}
}
}
#ifdef UNEQUALSOFTENINGS
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
h = All.ForceSofteningQ;
maxsofttype = (nop->u.d.bitflags >> 2) & 7;
if(maxsofttype == 7) /* may only occur for zero mass top-level nodes */
{
if(mass > 0)
endrun(986);
no = nop->u.d.nextnode;
continue;
}
else
{
if(h < All.ForceSoftening[maxsofttype])
{
h = All.ForceSoftening[maxsofttype];
if(r2 < h * h)
{
if(((nop->u.d.bitflags >> 5) & 1)) /* bit-5 signals that there are particles of different softening in the node */
{
no = nop->u.d.nextnode;
continue;
}
}
}
}
#else
if(ptype == 0)
h = soft;
else
h = All.ForceSofteningQ;
if(h < nop->maxsoft)
{
h = nop->maxsoft;
if(r2 < h * h)
{
no = nop->u.d.nextnode;
continue;
}
}
#endif
#endif
no = nop->u.d.sibling; /* ok, node can be used */
if(mode == 1)
{
if(((nop->u.d.bitflags) & 1)) /* Bit 0 signals that this node belongs to top-level tree */
continue;
}
}
r = sqrt(r2);
if(r >= h)
fac = mass / (r2 * r);
else
{
#ifdef UNEQUALSOFTENINGS
h_inv = 1.0 / h;
h3_inv = h_inv * h_inv * h_inv;
#endif
u = r * h_inv;
if(u < 0.5)
fac = mass * h3_inv * (10.666666666667 + u * u * (32.0 * u - 38.4));
else
fac =
mass * h3_inv * (21.333333333333 - 48.0 * u +
38.4 * u * u - 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
}
acc_x += dx * fac;
acc_y += dy * fac;
acc_z += dz * fac;
ninteractions++;
}
/* store result at the proper place */
if(mode == 0)
{
Q[target].GravAccel[0] = acc_x;
Q[target].GravAccel[1] = acc_y;
Q[target].GravAccel[2] = acc_z;
Q[target].GravCost = ninteractions;
}
else
{
GravDataResult[target].u.Acc[0] = acc_x;
GravDataResult[target].u.Acc[1] = acc_y;
GravDataResult[target].u.Acc[2] = acc_z;
GravDataResult[target].w.Ninteractions = ninteractions;
}
#ifdef PERIODIC
*ewaldcountsum += force_treeevaluate_ewald_correction(target, mode, pos_x, pos_y, pos_z, aold);
#endif
return ninteractions;
}
#ifdef PMGRID
/*! In the TreePM algorithm, the tree is walked only locally around the
* target coordinate. Tree nodes that fall outside a box of half
* side-length Rcut= RCUT*ASMTH*MeshSize can be discarded. The short-range
* potential is modified by a complementary error function, multiplied
* with the Newtonian form. The resulting short-range suppression compared
* to the Newtonian force is tabulated, because looking up from this table
* is faster than recomputing the corresponding factor, despite the
* memory-access panelty (which reduces cache performance) incurred by the
* table.
*/
int force_treeevaluate_shortrange_sub(int target, int mode)
{
struct NODE *nop = 0;
int no, ptype, ninteractions, tabindex;
double r2, dx, dy, dz, mass, r, fac, u, h, h_inv, h3_inv;
double acc_x, acc_y, acc_z, pos_x, pos_y, pos_z, aold;
double eff_dist;
double rcut, asmth, asmthfac, rcut2, dist;
#if defined(UNEQUALSOFTENINGS) && !defined(ADAPTIVE_GRAVSOFT_FORGAS)
int maxsofttype;
#endif
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
double soft = 0;
#endif
#ifdef PERIODIC
double boxsize, boxhalf;
boxsize = All.BoxSize;
boxhalf = 0.5 * All.BoxSize;
#endif
acc_x = 0;
acc_y = 0;
acc_z = 0;
ninteractions = 0;
if(mode == 0)
{
pos_x = Q[target].Pos[0];
pos_y = Q[target].Pos[1];
pos_z = Q[target].Pos[2];
ptype = Q[target].Type;
aold = All.ErrTolForceAcc * Q[target].OldAcc;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
soft = SphQ[target].Hsml;
#endif
}
else
{
pos_x = GravDataGet[target].u.Pos[0];
pos_y = GravDataGet[target].u.Pos[1];
pos_z = GravDataGet[target].u.Pos[2];
#ifdef UNEQUALSOFTENINGS
ptype = GravDataGet[target].Type;
#else
ptype = P[0].Type;
#endif
aold = All.ErrTolForceAcc * GravDataGet[target].w.OldAcc;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
soft = GravDataGet[target].Soft;
#endif
}
rcut = All.Rcut[0];
asmth = All.Asmth[0];
#ifdef PLACEHIGHRESREGION
if(((1 << ptype) & (PLACEHIGHRESREGION)))
{
rcut = All.Rcut[1];
asmth = All.Asmth[1];
}
#endif
rcut2 = rcut * rcut;
asmthfac = 0.5 / asmth * (NTAB / 3.0);
#ifndef UNEQUALSOFTENINGS
h = All.ForceSofteningQ;
h_inv = 1.0 / h;
h3_inv = h_inv * h_inv * h_inv;
#endif
no = All.MaxPart; /* root node */
while(no >= 0)
{
if(no < All.MaxPart)
{
/* the index of the node is the index of the particle */
dx = P[no].Pos[0] - pos_x;
dy = P[no].Pos[1] - pos_y;
dz = P[no].Pos[2] - pos_z;
#ifdef PERIODIC
dx = NEAREST(dx);
dy = NEAREST(dy);
dz = NEAREST(dz);
#endif
r2 = dx * dx + dy * dy + dz * dz;
mass = P[no].Mass;
#ifdef UNEQUALSOFTENINGS
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
h = soft;
else
h = All.ForceSofteningQ;
if(P[no].Type == 0)
{
if(h < SphP[no].Hsml)
h = SphP[no].Hsml;
}
else
{
if(h < All.ForceSoftening[P[no].Type])
h = All.ForceSoftening[P[no].Type];
}
#else
h = All.ForceSofteningQ;
if(h < All.ForceSoftening[P[no].Type])
h = All.ForceSoftening[P[no].Type];
#endif
#endif
no = Nextnode[no];
}
else /* we have an internal node */
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
if(mode == 0)
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
}
no = Nextnode[no - MaxNodes];
continue;
}
nop = &Nodes[no];
if(mode == 1)
{
if((nop->u.d.bitflags & 3) == 1) /* if it's a top-level node
* which does not contain
* local particles we can
* continue at this point
*/
{
no = nop->u.d.sibling;
continue;
}
}
mass = nop->u.d.mass;
dx = nop->u.d.s[0] - pos_x;
dy = nop->u.d.s[1] - pos_y;
dz = nop->u.d.s[2] - pos_z;
#ifdef PERIODIC
dx = NEAREST(dx);
dy = NEAREST(dy);
dz = NEAREST(dz);
#endif
r2 = dx * dx + dy * dy + dz * dz;
if(r2 > rcut2)
{
/* check whether we can stop walking along this branch */
eff_dist = rcut + 0.5 * nop->len;
#ifdef PERIODIC
dist = NEAREST(nop->center[0] - pos_x);
#else
dist = nop->center[0] - pos_x;
#endif
if(dist < -eff_dist || dist > eff_dist)
{
no = nop->u.d.sibling;
continue;
}
#ifdef PERIODIC
dist = NEAREST(nop->center[1] - pos_y);
#else
dist = nop->center[1] - pos_y;
#endif
if(dist < -eff_dist || dist > eff_dist)
{
no = nop->u.d.sibling;
continue;
}
#ifdef PERIODIC
dist = NEAREST(nop->center[2] - pos_z);
#else
dist = nop->center[2] - pos_z;
#endif
if(dist < -eff_dist || dist > eff_dist)
{
no = nop->u.d.sibling;
continue;
}
}
if(All.ErrTolTheta) /* check Barnes-Hut opening criterion */
{
if(nop->len * nop->len > r2 * All.ErrTolTheta * All.ErrTolTheta)
{
/* open cell */
no = nop->u.d.nextnode;
continue;
}
}
else /* check relative opening criterion */
{
if(mass * nop->len * nop->len > r2 * r2 * aold)
{
/* open cell */
no = nop->u.d.nextnode;
continue;
}
/* check in addition whether we lie inside the cell */
if(fabs(nop->center[0] - pos_x) < 0.60 * nop->len)
{
if(fabs(nop->center[1] - pos_y) < 0.60 * nop->len)
{
if(fabs(nop->center[2] - pos_z) < 0.60 * nop->len)
{
no = nop->u.d.nextnode;
continue;
}
}
}
}
#ifdef UNEQUALSOFTENINGS
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
h = All.ForceSofteningQ;
maxsofttype = (nop->u.d.bitflags >> 2) & 7;
if(maxsofttype == 7) /* may only occur for zero mass top-level nodes */
{
if(mass > 0)
endrun(987);
no = nop->u.d.nextnode;
continue;
}
else
{
if(h < All.ForceSoftening[maxsofttype])
{
h = All.ForceSoftening[maxsofttype];
if(r2 < h * h)
{
if(((nop->u.d.bitflags >> 5) & 1)) /* bit-5 signals that there are particles of different softening in the node */
{
no = nop->u.d.nextnode;
continue;
}
}
}
}
#else
if(ptype == 0)
h = soft;
else
h = All.ForceSofteningQ;
if(h < nop->maxsoft)
{
h = nop->maxsoft;
if(r2 < h * h)
{
no = nop->u.d.nextnode;
continue;
}
}
#endif
#endif
no = nop->u.d.sibling; /* ok, node can be used */
if(mode == 1)
{
if(((nop->u.d.bitflags) & 1)) /* Bit 0 signals that this node belongs to top-level tree */
continue;
}
}
r = sqrt(r2);
if(r >= h)
fac = mass / (r2 * r);
else
{
#ifdef UNEQUALSOFTENINGS
h_inv = 1.0 / h;
h3_inv = h_inv * h_inv * h_inv;
#endif
u = r * h_inv;
if(u < 0.5)
fac = mass * h3_inv * (10.666666666667 + u * u * (32.0 * u - 38.4));
else
fac =
mass * h3_inv * (21.333333333333 - 48.0 * u +
38.4 * u * u - 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
}
tabindex = (int) (asmthfac * r);
if(tabindex < NTAB)
{
fac *= shortrange_table[tabindex];
acc_x += dx * fac;
acc_y += dy * fac;
acc_z += dz * fac;
ninteractions++;
}
}
/* store result at the proper place */
if(mode == 0)
{
Q[target].GravAccel[0] = acc_x;
Q[target].GravAccel[1] = acc_y;
Q[target].GravAccel[2] = acc_z;
Q[target].GravCost = ninteractions;
}
else
{
GravDataResult[target].u.Acc[0] = acc_x;
GravDataResult[target].u.Acc[1] = acc_y;
GravDataResult[target].u.Acc[2] = acc_z;
GravDataResult[target].w.Ninteractions = ninteractions;
}
return ninteractions;
}
#endif
/*! This routine computes the gravitational potential by walking the
* tree. The same opening criteria is used as for the gravitational force
* walk.
*/
void force_treeevaluate_potential_sub(int target, int mode)
{
struct NODE *nop = 0;
int no, ptype;
double r2, dx, dy, dz, mass, r, u, h, h_inv, wp;
double pot, pos_x, pos_y, pos_z, aold;
#if defined(UNEQUALSOFTENINGS) && !defined(ADAPTIVE_GRAVSOFT_FORGAS)
int maxsofttype;
#endif
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
double soft = 0;
#endif
#ifdef PERIODIC
double boxsize, boxhalf;
boxsize = All.BoxSize;
boxhalf = 0.5 * All.BoxSize;
#endif
pot = 0;
if(mode == 0)
{
pos_x = Q[target].Pos[0];
pos_y = Q[target].Pos[1];
pos_z = Q[target].Pos[2];
ptype = Q[target].Type;
aold = All.ErrTolForceAcc * Q[target].OldAcc;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
soft = SphQ[target].Hsml;
#endif
}
else
{
pos_x = GravDataGet[target].u.Pos[0];
pos_y = GravDataGet[target].u.Pos[1];
pos_z = GravDataGet[target].u.Pos[2];
#ifdef UNEQUALSOFTENINGS
ptype = GravDataGet[target].Type;
#else
ptype = P[0].Type;
#endif
aold = All.ErrTolForceAcc * GravDataGet[target].w.OldAcc;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
soft = GravDataGet[target].Soft;
#endif
}
#ifndef UNEQUALSOFTENINGS
h = All.ForceSofteningQ;
h_inv = 1.0 / h;
#endif
no = All.MaxPart;
while(no >= 0)
{
if(no < All.MaxPart) /* single particle */
{
/* the index of the node is the index of the particle */
/* observe the sign */
dx = P[no].Pos[0] - pos_x;
dy = P[no].Pos[1] - pos_y;
dz = P[no].Pos[2] - pos_z;
mass = P[no].Mass;
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
if(mode == 0)
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
}
no = Nextnode[no - MaxNodes];
continue;
}
nop = &Nodes[no];
dx = nop->u.d.s[0] - pos_x;
dy = nop->u.d.s[1] - pos_y;
dz = nop->u.d.s[2] - pos_z;
mass = nop->u.d.mass;
}
#ifdef PERIODIC
dx = NEAREST(dx);
dy = NEAREST(dy);
dz = NEAREST(dz);
#endif
r2 = dx * dx + dy * dy + dz * dz;
if(no < All.MaxPart)
{
#ifdef UNEQUALSOFTENINGS
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
h = soft;
else
h = All.ForceSofteningQ;
if(P[no].Type == 0)
{
if(h < SphP[no].Hsml)
h = SphP[no].Hsml;
}
else
{
if(h < All.ForceSoftening[P[no].Type])
h = All.ForceSoftening[P[no].Type];
}
#else
h = All.ForceSofteningQ;
if(h < All.ForceSoftening[P[no].Type])
h = All.ForceSoftening[P[no].Type];
#endif
#endif
no = Nextnode[no];
}
else /* we have an internal node. Need to check opening criterion */
{
if(mode == 1)
{
if((nop->u.d.bitflags & 3) == 1) /* if it's a top-level node
* which does not contain
* local particles we can make
* a short-cut
*/
{
no = nop->u.d.sibling;
continue;
}
}
if(All.ErrTolTheta) /* check Barnes-Hut opening criterion */
{
if(nop->len * nop->len > r2 * All.ErrTolTheta * All.ErrTolTheta)
{
/* open cell */
no = nop->u.d.nextnode;
continue;
}
}
else /* check relative opening criterion */
{
if(mass * nop->len * nop->len > r2 * r2 * aold)
{
/* open cell */
no = nop->u.d.nextnode;
continue;
}
if(fabs(nop->center[0] - pos_x) < 0.60 * nop->len)
{
if(fabs(nop->center[1] - pos_y) < 0.60 * nop->len)
{
if(fabs(nop->center[2] - pos_z) < 0.60 * nop->len)
{
no = nop->u.d.nextnode;
continue;
}
}
}
}
#ifdef UNEQUALSOFTENINGS
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
h = All.ForceSofteningQ;
maxsofttype = (nop->u.d.bitflags >> 2) & 7;
if(maxsofttype == 7) /* may only occur for zero mass top-level nodes */
{
if(mass > 0)
endrun(988);
no = nop->u.d.nextnode;
continue;
}
else
{
if(h < All.ForceSoftening[maxsofttype])
{
h = All.ForceSoftening[maxsofttype];
if(r2 < h * h)
{
if(((nop->u.d.bitflags >> 5) & 1)) /* bit-5 signals that there are particles of different softening in the node */
{
no = nop->u.d.nextnode;
continue;
}
}
}
}
#else
if(ptype == 0)
h = soft;
else
h = All.ForceSofteningQ;
if(h < nop->maxsoft)
{
h = nop->maxsoft;
if(r2 < h * h)
{
no = nop->u.d.nextnode;
continue;
}
}
#endif
#endif
no = nop->u.d.sibling; /* node can be used */
if(mode == 1)
{
if(((nop->u.d.bitflags) & 1)) /* Bit 0 signals that this node belongs to top-level tree */
continue;
}
}
r = sqrt(r2);
if(r >= h)
pot -= mass / r;
else
{
#ifdef UNEQUALSOFTENINGS
h_inv = 1.0 / h;
#endif
u = r * h_inv;
if(u < 0.5)
wp = -2.8 + u * u * (5.333333333333 + u * u * (6.4 * u - 9.6));
else
wp =
-3.2 + 0.066666666667 / u + u * u * (10.666666666667 +
u * (-16.0 + u * (9.6 - 2.133333333333 * u)));
pot += mass * h_inv * wp;
}
#ifdef PERIODIC
pot += mass * ewald_pot_corr(dx, dy, dz);
#endif
}
/* store result at the proper place */
if(mode == 0)
Q[target].Potential = pot;
else
GravDataResult[target].u.Potential = pot;
}
#ifdef PMGRID
/*! This function computes the short-range potential when the TreePM
* algorithm is used. This potential is the Newtonian potential, modified
* by a complementary error function.
*/
void force_treeevaluate_potential_shortrange_sub(int target, int mode)
{
struct NODE *nop = 0;
int no, ptype, tabindex;
double r2, dx, dy, dz, mass, r, u, h, h_inv, wp;
double pot, pos_x, pos_y, pos_z, aold;
double eff_dist, fac, rcut, asmth, asmthfac;
double dxx, dyy, dzz;
#if defined(UNEQUALSOFTENINGS) && !defined(ADAPTIVE_GRAVSOFT_FORGAS)
int maxsofttype;
#endif
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
double soft = 0;
#endif
#ifdef PERIODIC
double boxsize, boxhalf;
boxsize = All.BoxSize;
boxhalf = 0.5 * All.BoxSize;
#endif
pot = 0;
if(mode == 0)
{
pos_x = Q[target].Pos[0];
pos_y = Q[target].Pos[1];
pos_z = Q[target].Pos[2];
ptype = Q[target].Type;
aold = All.ErrTolForceAcc * Q[target].OldAcc;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
soft = SphQ[target].Hsml;
#endif
}
else
{
pos_x = GravDataGet[target].u.Pos[0];
pos_y = GravDataGet[target].u.Pos[1];
pos_z = GravDataGet[target].u.Pos[2];
#ifdef UNEQUALSOFTENINGS
ptype = GravDataGet[target].Type;
#else
ptype = P[0].Type;
#endif
aold = All.ErrTolForceAcc * GravDataGet[target].w.OldAcc;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
soft = GravDataGet[target].Soft;
#endif
}
rcut = All.Rcut[0];
asmth = All.Asmth[0];
#ifdef PLACEHIGHRESREGION
if(((1 << ptype) & (PLACEHIGHRESREGION)))
{
rcut = All.Rcut[1];
asmth = All.Asmth[1];
}
#endif
asmthfac = 0.5 / asmth * (NTAB / 3.0);
#ifndef UNEQUALSOFTENINGS
h = All.ForceSofteningQ;
h_inv = 1.0 / h;
#endif
no = All.MaxPart;
while(no >= 0)
{
if(no < All.MaxPart) /* single particle */
{
/* the index of the node is the index of the particle */
/* observe the sign */
dx = P[no].Pos[0] - pos_x;
dy = P[no].Pos[1] - pos_y;
dz = P[no].Pos[2] - pos_z;
mass = P[no].Mass;
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
if(mode == 0)
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
}
no = Nextnode[no - MaxNodes];
continue;
}
nop = &Nodes[no];
dx = nop->u.d.s[0] - pos_x;
dy = nop->u.d.s[1] - pos_y;
dz = nop->u.d.s[2] - pos_z;
mass = nop->u.d.mass;
}
#ifdef PERIODIC
dx = NEAREST(dx);
dy = NEAREST(dy);
dz = NEAREST(dz);
#endif
r2 = dx * dx + dy * dy + dz * dz;
if(no < All.MaxPart)
{
#ifdef UNEQUALSOFTENINGS
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
h = soft;
else
h = All.ForceSofteningQ;
if(P[no].Type == 0)
{
if(h < SphP[no].Hsml)
h = SphP[no].Hsml;
}
else
{
if(h < All.ForceSoftening[P[no].Type])
h = All.ForceSoftening[P[no].Type];
}
#else
h = All.ForceSofteningQ;
if(h < All.ForceSoftening[P[no].Type])
h = All.ForceSoftening[P[no].Type];
#endif
#endif
no = Nextnode[no];
}
else /* we have an internal node. Need to check opening criterion */
{
/* check whether we can stop walking along this branch */
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
if(mode == 0)
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
}
no = Nextnode[no - MaxNodes];
continue;
}
if(mode == 1)
{
if((nop->u.d.bitflags & 3) == 1) /* if it's a top-level node which does not contain local particles */
{
no = nop->u.d.sibling;
continue;
}
}
eff_dist = rcut + 0.5 * nop->len;
dxx = nop->center[0] - pos_x; /* observe the sign ! */
dyy = nop->center[1] - pos_y; /* this vector is -y in my thesis notation */
dzz = nop->center[2] - pos_z;
#ifdef PERIODIC
dxx = NEAREST(dxx);
dyy = NEAREST(dyy);
dzz = NEAREST(dzz);
#endif
if(dxx < -eff_dist || dxx > eff_dist)
{
no = nop->u.d.sibling;
continue;
}
if(dyy < -eff_dist || dyy > eff_dist)
{
no = nop->u.d.sibling;
continue;
}
if(dzz < -eff_dist || dzz > eff_dist)
{
no = nop->u.d.sibling;
continue;
}
if(All.ErrTolTheta) /* check Barnes-Hut opening criterion */
{
if(nop->len * nop->len > r2 * All.ErrTolTheta * All.ErrTolTheta)
{
/* open cell */
no = nop->u.d.nextnode;
continue;
}
}
else /* check relative opening criterion */
{
if(mass * nop->len * nop->len > r2 * r2 * aold)
{
/* open cell */
no = nop->u.d.nextnode;
continue;
}
if(fabs(nop->center[0] - pos_x) < 0.60 * nop->len)
{
if(fabs(nop->center[1] - pos_y) < 0.60 * nop->len)
{
if(fabs(nop->center[2] - pos_z) < 0.60 * nop->len)
{
no = nop->u.d.nextnode;
continue;
}
}
}
}
#ifdef UNEQUALSOFTENINGS
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
h = All.ForceSofteningQ;
maxsofttype = (nop->u.d.bitflags >> 2) & 7;
if(maxsofttype == 7) /* may only occur for zero mass top-level nodes */
{
if(mass > 0)
endrun(989);
no = nop->u.d.nextnode;
continue;
}
else
{
if(h < All.ForceSoftening[maxsofttype])
{
h = All.ForceSoftening[maxsofttype];
if(r2 < h * h)
{
/* bit-5 signals that there are particles of
* different softening in the node
*/
if(((nop->u.d.bitflags >> 5) & 1))
{
no = nop->u.d.nextnode;
continue;
}
}
}
}
#else
if(ptype == 0)
h = soft;
else
h = All.ForceSofteningQ;
if(h < nop->maxsoft)
{
h = nop->maxsoft;
if(r2 < h * h)
{
no = nop->u.d.nextnode;
continue;
}
}
#endif
#endif
no = nop->u.d.sibling; /* node can be used */
if(mode == 1)
{
if(((nop->u.d.bitflags) & 1)) /* Bit 0 signals that this node belongs to top-level tree */
continue;
}
}
r = sqrt(r2);
tabindex = (int) (r * asmthfac);
if(tabindex < NTAB)
{
fac = shortrange_table_potential[tabindex];
if(r >= h)
pot -= fac * mass / r;
else
{
#ifdef UNEQUALSOFTENINGS
h_inv = 1.0 / h;
#endif
u = r * h_inv;
if(u < 0.5)
wp = -2.8 + u * u * (5.333333333333 + u * u * (6.4 * u - 9.6));
else
wp =
-3.2 + 0.066666666667 / u + u * u * (10.666666666667 +
u * (-16.0 + u * (9.6 - 2.133333333333 * u)));
pot += fac * mass * h_inv * wp;
}
}
}
/* store result at the proper place */
if(mode == 0)
Q[target].Potential = pot;
else
GravDataResult[target].u.Potential = pot;
}
#endif
#endif /*PY_INTERFACE*/
c4science · Help