Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85231815
forcetree_sub.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Sep 27, 16:36
Size
28 KB
Mime Type
text/x-c
Expires
Sun, Sep 29, 16:36 (2 d)
Engine
blob
Format
Raw Data
Handle
21145928
Attached To
rGEAR Gear
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*/
Event Timeline
Log In to Comment