Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F106897661
sticky.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
Wed, Apr 2, 05:52
Size
30 KB
Mime Type
text/x-c
Expires
Fri, Apr 4, 05:52 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
25300827
Attached To
rPNBODY pNbody
sticky.c
View Options
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include "allvars.h"
#include "proto.h"
#ifdef PERIODIC
static
double
boxSize
,
boxHalf
;
#ifdef LONG_X
static
double
boxSize_X
,
boxHalf_X
;
#else
#define boxSize_X boxSize
#define boxHalf_X boxHalf
#endif
#ifdef LONG_Y
static
double
boxSize_Y
,
boxHalf_Y
;
#else
#define boxSize_Y boxSize
#define boxHalf_Y boxHalf
#endif
#ifdef LONG_Z
static
double
boxSize_Z
,
boxHalf_Z
;
#else
#define boxSize_Z boxSize
#define boxHalf_Z boxHalf
#endif
#endif
#ifdef MULTIPHASE
static
int
Ncz
;
static
int
Ncxy
;
static
int
Ncxyz
;
static
double
hubble_a
,
atime
,
hubble_a2
,
fac_mu
,
fac_vsic_fix
,
a3inv
,
fac_egy
;
/*! \file sticky.c
* \brief Compute sticky collisions
*
*/
/*! init sticky
*
*/
void
init_sticky
()
{
Ncz
=
All
.
StickyGridNz
;
Ncxy
=
All
.
StickyGridNx
*
All
.
StickyGridNy
;
Ncxyz
=
All
.
StickyGridNx
*
All
.
StickyGridNy
*
All
.
StickyGridNz
;
}
/*! This function is the driver routine for the calculation of sticky collisions
*/
void
sticky
()
{
double
t0
,
t1
;
double
dt
;
int
Tis
,
Tic
;
t0
=
second
();
/* measure the time for the full chimie computation */
Tis
=
log
(
All
.
StickyLastCollisionTime
/
All
.
TimeBegin
)
/
All
.
Timebase_interval
;
Tic
=
All
.
Ti_Current
;
dt
=
get_cosmictime_difference
(
Tis
,
Tic
);
//if((All.Time - All.StickyLastCollisionTime) >= All.StickyCollisionTime)
if
(
dt
>=
All
.
StickyCollisionTime
)
{
sticky_compute_energy_kin
(
1
);
//sticky_collisions(); /* old implementation */
if
(
ThisTask
==
0
)
printf
(
"sticky first loop
\n
"
);
sticky_collisions2
(
1
);
if
(
ThisTask
==
0
)
printf
(
"sticky second loop
\n
"
);
sticky_collisions2
(
2
);
sticky_compute_energy_kin
(
2
);
All
.
StickyLastCollisionTime
+=
All
.
StickyCollisionTime
;
}
t1
=
second
();
All
.
CPU_Sticky
+=
timediff
(
t0
,
t1
);
}
/*! This function is used as a comparison kernel in a sort routine.
*/
int
compare_sticky_index
(
const
void
*
a
,
const
void
*
b
)
{
if
(((
struct
Sticky_index
*
)
a
)
->
CellIndex
<
(((
struct
Sticky_index
*
)
b
)
->
CellIndex
))
return
-
1
;
if
(((
struct
Sticky_index
*
)
a
)
->
CellIndex
>
(((
struct
Sticky_index
*
)
b
)
->
CellIndex
))
return
+
1
;
return
0
;
}
/*! This function returns the first particle that may sticky-interact with
* the given particle. This routine use a pseudo-grid to find pairs, based
* on the Bournaud technique
*/
#ifdef MULTIPHASE
int
find_sticky_pairs
(
int
ii
)
{
int
jj
;
int
ci
,
cj
;
int
i
,
j
;
double
dx
,
dy
,
dz
,
r2
,
r
;
double
dvx
,
dvy
,
dvz
,
vdotr
,
vdotr2
,
dv
;
if
(
ii
>=
N_gas
)
return
-
11
;
ci
=
StickyIndex
[
ii
].
CellIndex
;
if
(
ci
>=
Ncxyz
)
return
-
12
;
jj
=
ii
;
do
{
jj
++
;
if
(
jj
>=
N_gas
)
return
-
1
;
cj
=
StickyIndex
[
jj
].
CellIndex
;
if
(
cj
>=
Ncxyz
)
return
-
2
;
if
(
ci
!=
cj
)
return
-
3
;
if
(
SphP
[
StickyIndex
[
jj
].
Index
].
StickyFlag
)
/* if its free, return it, else, go to the next one */
{
/* here we check that we can interact */
i
=
StickyIndex
[
ii
].
Index
;
j
=
StickyIndex
[
jj
].
Index
;
dx
=
P
[
i
].
Pos
[
0
]
-
P
[
j
].
Pos
[
0
];
dy
=
P
[
i
].
Pos
[
1
]
-
P
[
j
].
Pos
[
1
];
dz
=
P
[
i
].
Pos
[
2
]
-
P
[
j
].
Pos
[
2
];
r2
=
dx
*
dx
+
dy
*
dy
+
dz
*
dz
;
r
=
sqrt
(
r2
);
dvx
=
P
[
i
].
Vel
[
0
]
-
P
[
j
].
Vel
[
0
];
dvy
=
P
[
i
].
Vel
[
1
]
-
P
[
j
].
Vel
[
1
];
dvz
=
P
[
i
].
Vel
[
2
]
-
P
[
j
].
Vel
[
2
];
vdotr
=
dx
*
dvx
+
dy
*
dvy
+
dz
*
dvz
;
vdotr2
=
vdotr
;
dv
=
vdotr2
/
r
;
if
(
(
vdotr2
<
0
)
&&
(
-
vdotr2
>
All
.
StickyMinVelocity
))
{
printf
(
"(1)%d %d
\n
"
,
i
,
j
);
return
StickyIndex
[
jj
].
Index
;
/* ok, we accept the collision */
}
}
}
while
(
1
);
}
#endif
void
sticky_compute_energy_kin
(
int
mode
)
{
int
i
;
double
DeltaEgyKin
;
double
Tot_DeltaEgyKin
;
DeltaEgyKin
=
0
;
Tot_DeltaEgyKin
=
0
;
if
(
mode
==
1
)
{
LocalSysState
.
EnergyKin1
=
0
;
LocalSysState
.
EnergyKin2
=
0
;
}
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
{
if
(
P
[
i
].
Type
==
0
)
{
if
(
mode
==
1
)
LocalSysState
.
EnergyKin1
+=
0.5
*
P
[
i
].
Mass
*
(
P
[
i
].
Vel
[
0
]
*
P
[
i
].
Vel
[
0
]
+
P
[
i
].
Vel
[
1
]
*
P
[
i
].
Vel
[
1
]
+
P
[
i
].
Vel
[
2
]
*
P
[
i
].
Vel
[
2
]);
else
LocalSysState
.
EnergyKin2
+=
0.5
*
P
[
i
].
Mass
*
(
P
[
i
].
Vel
[
0
]
*
P
[
i
].
Vel
[
0
]
+
P
[
i
].
Vel
[
1
]
*
P
[
i
].
Vel
[
1
]
+
P
[
i
].
Vel
[
2
]
*
P
[
i
].
Vel
[
2
]);
}
}
if
(
mode
==
2
)
{
DeltaEgyKin
=
LocalSysState
.
EnergyKin2
-
LocalSysState
.
EnergyKin1
;
MPI_Reduce
(
&
DeltaEgyKin
,
&
Tot_DeltaEgyKin
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
0
,
MPI_COMM_WORLD
);
LocalSysState
.
EnergyRadSticky
-=
DeltaEgyKin
;
}
}
/*! compute sticky collisions
*
*/
void
sticky_collisions
()
{
int
i
,
j
,
k
;
int
startnode
;
FLOAT
*
pos
,
*
vel
;
FLOAT
mass
,
h_i
;
int
phase
;
double
vcm
[
3
],
newv1
[
3
],
newv2
[
3
];
double
dv1
[
3
],
dv2
[
3
];
double
dx
,
dy
,
dz
,
dvx
,
dvy
,
dvz
;
double
r
,
r2
;
double
vdotr
,
vdotr2
;
double
dv
;
double
beta_r
,
beta_t
;
double
M1M
,
M2M
,
M12
;
double
dbeta_dv
,
dbeta_dv_er_x
,
dbeta_dv_er_y
,
dbeta_dv_er_z
;
double
dv_beta_t_x
,
dv_beta_t_y
,
dv_beta_t_z
;
NumColLocal
=
0
;
NumNoColLocal
=
0
;
beta_r
=
All
.
StickyBetaR
;
beta_t
=
All
.
StickyBetaT
;
double
xi
,
yi
,
zi
;
int
ix
,
iy
,
iz
;
size_t
bytes
;
if
(
ThisTask
==
0
)
printf
(
"Start sticky collisions...
\n
"
);
/*****************************/
/* first, update StickyIndex */
/*****************************/
if
(
All
.
StickyUseGridForCollisions
)
{
/* allocate memory */
if
(
!
(
StickyIndex
=
malloc
(
bytes
=
N_gas
*
sizeof
(
struct
Sticky_index
))))
{
printf
(
"failed to allocate memory for `StickyIndex' (%g MB).
\n
"
,
bytes
/
(
1024.0
*
1024.0
));
endrun
(
2233
);
}
/* loop gas particles */
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
{
//if (P[i].Type==0) /* we also take also into account the stars not transfered in stars now */
{
StickyIndex
[
i
].
Index
=
i
;
StickyIndex
[
i
].
CellIndex
=
Ncxyz
;
StickyIndex
[
i
].
Flag
=
1
;
if
(
SphP
[
i
].
StickyFlag
)
/* particle that may collide (set in phase.c) */
{
/* scale between 0 and nx,ny,nz */
xi
=
All
.
StickyGridNx
*
(
P
[
i
].
Pos
[
0
]
-
All
.
StickyGridXmin
)
/
(
All
.
StickyGridXmax
-
All
.
StickyGridXmin
);
yi
=
All
.
StickyGridNy
*
(
P
[
i
].
Pos
[
1
]
-
All
.
StickyGridYmin
)
/
(
All
.
StickyGridYmax
-
All
.
StickyGridYmin
);
zi
=
All
.
StickyGridNz
*
(
P
[
i
].
Pos
[
2
]
-
All
.
StickyGridZmin
)
/
(
All
.
StickyGridZmax
-
All
.
StickyGridZmin
);
ix
=
(
int
)(
xi
);
iy
=
(
int
)(
yi
);
iz
=
(
int
)(
zi
);
if
(
ix
>=
0
&&
ix
<
All
.
StickyGridNx
)
if
(
iy
>=
0
&&
iy
<
All
.
StickyGridNy
)
if
(
iz
>=
0
&&
iz
<
All
.
StickyGridNz
)
{
j
=
(
ix
)
*
(
Ncxy
)
+
(
iy
)
*
(
Ncz
)
+
(
iz
);
StickyIndex
[
i
].
Index
=
i
;
StickyIndex
[
i
].
CellIndex
=
j
;
}
}
}
}
/* sort particles */
qsort
(
StickyIndex
,
N_gas
,
sizeof
(
struct
Sticky_index
),
compare_sticky_index
);
/* record index in SphP (not necessary, but simplifies) */
/* loop over cells */
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
{
//printf("%d %d\n",StickyIndex[i].CellIndex,StickyIndex[i].Index);
j
=
StickyIndex
[
i
].
Index
;
SphP
[
j
].
StickyIndex
=
i
;
}
}
/*****************************/
/* loop over all particles */
/*****************************/
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
{
//#ifdef SFR
// if((P[i].Ti_endstep == All.Ti_Current) && (P[i].Type == 0))
//#else
// if(P[i].Ti_endstep == All.Ti_Current)
//#endif
#ifdef SFR
if
(
P
[
i
].
Type
==
0
)
#endif
{
if
(
SphP
[
i
].
Phase
==
GAS_STICKY
)
{
if
(
SphP
[
i
].
StickyFlag
)
{
pos
=
P
[
i
].
Pos
;
vel
=
SphP
[
i
].
VelPred
;
mass
=
P
[
i
].
Mass
;
h_i
=
SphP
[
i
].
Hsml
;
phase
=
SphP
[
i
].
Phase
;
startnode
=
All
.
MaxPart
;
if
(
All
.
StickyUseGridForCollisions
)
j
=
find_sticky_pairs
(
SphP
[
i
].
StickyIndex
);
else
j
=
ngb_treefind_sticky_collisions
(
&
pos
[
0
],
h_i
,
phase
,
&
startnode
);
/* check that particles are in the same cell */
/*
if (j>=0 && j!=i)
{
int i1,i2;
xi = All.StickyGridNx*(P[i].Pos[0]-All.StickyGridXmin)/(All.StickyGridXmax-All.StickyGridXmin);
yi = All.StickyGridNy*(P[i].Pos[1]-All.StickyGridYmin)/(All.StickyGridYmax-All.StickyGridYmin);
zi = All.StickyGridNz*(P[i].Pos[2]-All.StickyGridZmin)/(All.StickyGridZmax-All.StickyGridZmin);
ix = (int)(xi);
iy = (int)(yi);
iz = (int)(zi);
i1 = (ix)*(Ncxy) + (iy)*(Ncz) + (iz);
xi = All.StickyGridNx*(P[i].Pos[0]-All.StickyGridXmin)/(All.StickyGridXmax-All.StickyGridXmin);
yi = All.StickyGridNy*(P[i].Pos[1]-All.StickyGridYmin)/(All.StickyGridYmax-All.StickyGridYmin);
zi = All.StickyGridNz*(P[i].Pos[2]-All.StickyGridZmin)/(All.StickyGridZmax-All.StickyGridZmin);
ix = (int)(xi);
iy = (int)(yi);
iz = (int)(zi);
i2 = (ix)*(Ncxy) + (iy)*(Ncz) + (iz);
if (i1!=i2)
{
printf("\n\nSomething odd here!\n");
endrun("999888777");
}
}
*/
if
(
j
>=
0
&&
j
!=
i
)
/* particles may collide */
{
if
(
j
>=
N_gas
)
{
printf
(
"sticky_collisions : j=%d,N_gas=%d j>N_gas !!!"
,
j
,
N_gas
);
exit
(
-
1
);
}
dx
=
pos
[
0
]
-
P
[
j
].
Pos
[
0
];
dy
=
pos
[
1
]
-
P
[
j
].
Pos
[
1
];
dz
=
pos
[
2
]
-
P
[
j
].
Pos
[
2
];
r2
=
dx
*
dx
+
dy
*
dy
+
dz
*
dz
;
r
=
sqrt
(
r2
);
dvx
=
vel
[
0
]
-
SphP
[
j
].
VelPred
[
0
];
dvy
=
vel
[
1
]
-
SphP
[
j
].
VelPred
[
1
];
dvz
=
vel
[
2
]
-
SphP
[
j
].
VelPred
[
2
];
vdotr
=
dx
*
dvx
+
dy
*
dvy
+
dz
*
dvz
;
vdotr2
=
vdotr
;
dv
=
vdotr2
/
r
;
if
(
vdotr2
<
0
)
/* particles approaches !!! */
{
/* compute new velocities */
M12
=
mass
+
P
[
j
].
Mass
;
M1M
=
mass
/
M12
;
M2M
=
P
[
j
].
Mass
/
M12
;
/* compute velocity of the mass center */
for
(
k
=
0
;
k
<
3
;
k
++
)
vcm
[
k
]
=
(
mass
*
vel
[
k
]
+
P
[
j
].
Mass
*
SphP
[
j
].
VelPred
[
k
])
/
M12
;
dbeta_dv
=
(
beta_r
-
beta_t
)
*
dv
;
dbeta_dv_er_x
=
-
dbeta_dv
*
dx
/
r
;
dbeta_dv_er_y
=
-
dbeta_dv
*
dy
/
r
;
dbeta_dv_er_z
=
-
dbeta_dv
*
dz
/
r
;
dv_beta_t_x
=
dvx
*
beta_t
;
dv_beta_t_y
=
dvy
*
beta_t
;
dv_beta_t_z
=
dvz
*
beta_t
;
newv1
[
0
]
=
M2M
*
(
+
dv_beta_t_x
-
dbeta_dv_er_x
)
+
vcm
[
0
];
newv1
[
1
]
=
M2M
*
(
+
dv_beta_t_y
-
dbeta_dv_er_y
)
+
vcm
[
1
];
newv1
[
2
]
=
M2M
*
(
+
dv_beta_t_z
-
dbeta_dv_er_z
)
+
vcm
[
2
];
newv2
[
0
]
=
M1M
*
(
-
dv_beta_t_x
+
dbeta_dv_er_x
)
+
vcm
[
0
];
newv2
[
1
]
=
M1M
*
(
-
dv_beta_t_y
+
dbeta_dv_er_y
)
+
vcm
[
1
];
newv2
[
2
]
=
M1M
*
(
-
dv_beta_t_z
+
dbeta_dv_er_z
)
+
vcm
[
2
];
/* new velocities */
for
(
k
=
0
;
k
<
3
;
k
++
)
{
dv1
[
k
]
=
newv1
[
k
]
-
SphP
[
i
].
VelPred
[
k
];
dv2
[
k
]
=
newv2
[
k
]
-
SphP
[
j
].
VelPred
[
k
];
SphP
[
i
].
VelPred
[
k
]
=
newv1
[
k
];
SphP
[
j
].
VelPred
[
k
]
=
newv2
[
k
];
P
[
i
].
Vel
[
k
]
+=
dv1
[
k
];
P
[
j
].
Vel
[
k
]
+=
dv2
[
k
];
}
/* set particles as non sticky-active */
SphP
[
i
].
StickyFlag
=
0
;
SphP
[
j
].
StickyFlag
=
0
;
SphP
[
i
].
StickyTime
=
All
.
Time
+
All
.
StickyIdleTime
;
SphP
[
j
].
StickyTime
=
All
.
Time
+
All
.
StickyIdleTime
;
/* record collision */
NumColLocal
+=
2
;
}
}
}
}
}
}
/* write some statistics */
MPI_Allreduce
(
&
NumColLocal
,
&
NumCol
,
1
,
MPI_INT
,
MPI_SUM
,
MPI_COMM_WORLD
);
MPI_Allreduce
(
&
NumColPotLocal
,
&
NumColPot
,
1
,
MPI_INT
,
MPI_SUM
,
MPI_COMM_WORLD
);
MPI_Allreduce
(
&
NumNoColLocal
,
&
NumNoCol
,
1
,
MPI_INT
,
MPI_SUM
,
MPI_COMM_WORLD
);
if
(
ThisTask
==
0
)
{
fprintf
(
FdSticky
,
"Step %d, Time: %g NumColPot: %d NumCol: %d NumNoCol: %d
\n
"
,
All
.
NumCurrentTiStep
,
All
.
Time
,(
int
)
NumColPot
,(
int
)
NumCol
,(
int
)
NumNoCol
);
fflush
(
FdSticky
);
}
if
(
All
.
StickyUseGridForCollisions
)
free
(
StickyIndex
);
if
(
ThisTask
==
0
)
printf
(
"sticky collisions done.
\n
"
);
}
/*! This function is the driver routine for the calculation of hydrodynamical
* force and rate of change of entropy due to shock heating for all active
* particles .
*
* During the first loop, each particle compute j (a potentially colliding particle)
* During the second loop, it check if it can collide and collide if it can.
*
*/
void
sticky_collisions2
(
int
loop
)
{
long
long
ntot
,
ntotleft
;
int
i
,
j
,
k
,
n
,
ngrp
,
maxfill
,
source
,
ndone
;
int
*
nbuffer
,
*
noffset
,
*
nsend_local
,
*
nsend
,
*
numlist
,
*
ndonelist
;
int
level
,
sendTask
,
recvTask
,
nexport
,
place
;
double
tstart
,
tend
,
sumt
,
sumcomm
;
double
timecomp
=
0
,
timecommsumm
=
0
,
timeimbalance
=
0
,
sumimbalance
;
MPI_Status
status
;
double
dv
[
3
];
#ifdef PERIODIC
boxSize
=
All
.
BoxSize
;
boxHalf
=
0.5
*
All
.
BoxSize
;
#ifdef LONG_X
boxHalf_X
=
boxHalf
*
LONG_X
;
boxSize_X
=
boxSize
*
LONG_X
;
#endif
#ifdef LONG_Y
boxHalf_Y
=
boxHalf
*
LONG_Y
;
boxSize_Y
=
boxSize
*
LONG_Y
;
#endif
#ifdef LONG_Z
boxHalf_Z
=
boxHalf
*
LONG_Z
;
boxSize_Z
=
boxSize
*
LONG_Z
;
#endif
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
double
v1m
,
v2m
;
#endif
if
(
All
.
ComovingIntegrationOn
)
{
/* Factors for comoving integration of hydro */
hubble_a
=
All
.
Omega0
/
(
All
.
Time
*
All
.
Time
*
All
.
Time
)
+
(
1
-
All
.
Omega0
-
All
.
OmegaLambda
)
/
(
All
.
Time
*
All
.
Time
)
+
All
.
OmegaLambda
;
hubble_a
=
All
.
Hubble
*
sqrt
(
hubble_a
);
hubble_a2
=
All
.
Time
*
All
.
Time
*
hubble_a
;
fac_mu
=
pow
(
All
.
Time
,
3
*
(
GAMMA
-
1
)
/
2
)
/
All
.
Time
;
fac_egy
=
pow
(
All
.
Time
,
3
*
(
GAMMA
-
1
));
fac_vsic_fix
=
hubble_a
*
pow
(
All
.
Time
,
3
*
GAMMA_MINUS1
);
a3inv
=
1
/
(
All
.
Time
*
All
.
Time
*
All
.
Time
);
atime
=
All
.
Time
;
#ifdef FEEDBACK
fac_pow
=
fac_egy
*
atime
*
atime
;
#endif
}
else
{
hubble_a
=
hubble_a2
=
atime
=
fac_mu
=
fac_vsic_fix
=
a3inv
=
fac_egy
=
1.0
;
#ifdef FEEDBACK
fac_pow
=
1.0
;
#endif
}
/* `NumSphUpdate' gives the number of particles on this processor that want a force update */
for
(
n
=
0
,
NumSphUpdate
=
0
;
n
<
N_gas
;
n
++
)
{
#ifdef SFR
//if((P[n].Ti_endstep == All.Ti_Current) && (P[n].Type == 0))
#else
//if(P[n].Ti_endstep == All.Ti_Current)
#endif
if
((
SphP
[
n
].
Phase
==
GAS_STICKY
)
&&
(
SphP
[
n
].
StickyFlag
)
)
{
if
(
loop
==
1
)
{
SphP
[
n
].
StickyMaxFs
=
0
;
SphP
[
n
].
StickyMaxID
=-
1
;
}
SphP
[
n
].
StickyNgb
=-
1
;
NumSphUpdate
++
;
}
}
numlist
=
malloc
(
NTask
*
sizeof
(
int
)
*
NTask
);
MPI_Allgather
(
&
NumSphUpdate
,
1
,
MPI_INT
,
numlist
,
1
,
MPI_INT
,
MPI_COMM_WORLD
);
for
(
i
=
0
,
ntot
=
0
;
i
<
NTask
;
i
++
)
ntot
+=
numlist
[
i
];
free
(
numlist
);
noffset
=
malloc
(
sizeof
(
int
)
*
NTask
);
/* offsets of bunches in common list */
nbuffer
=
malloc
(
sizeof
(
int
)
*
NTask
);
nsend_local
=
malloc
(
sizeof
(
int
)
*
NTask
);
nsend
=
malloc
(
sizeof
(
int
)
*
NTask
*
NTask
);
ndonelist
=
malloc
(
sizeof
(
int
)
*
NTask
);
i
=
0
;
/* first particle for this task */
ntotleft
=
ntot
;
/* particles left for all tasks together */
NumColLocal
=
0
;
NumNoColLocal
=
0
;
while
(
ntotleft
>
0
)
{
for
(
j
=
0
;
j
<
NTask
;
j
++
)
nsend_local
[
j
]
=
0
;
/* do local particles and prepare export list */
tstart
=
second
();
for
(
nexport
=
0
,
ndone
=
0
;
i
<
N_gas
&&
nexport
<
All
.
BunchSizeSticky
-
NTask
;
i
++
)
{
#ifdef SFR
//if((P[i].Ti_endstep == All.Ti_Current) && (P[i].Type == 0))
#else
//if(P[i].Ti_endstep == All.Ti_Current)
#endif
if
((
SphP
[
i
].
Phase
==
GAS_STICKY
)
&&
(
SphP
[
i
].
StickyFlag
))
{
ndone
++
;
for
(
j
=
0
;
j
<
NTask
;
j
++
)
Exportflag
[
j
]
=
0
;
sticky_evaluate
(
i
,
0
,
loop
);
/* the particle is not exported if an ngb has been found locally */
if
(
SphP
[
i
].
StickyNgb
!=-
1
)
for
(
j
=
0
;
j
<
NTask
;
j
++
)
Exportflag
[
j
]
=
0
;
for
(
j
=
0
;
j
<
NTask
;
j
++
)
{
if
(
Exportflag
[
j
])
{
for
(
k
=
0
;
k
<
3
;
k
++
)
{
StickyDataIn
[
nexport
].
Pos
[
k
]
=
P
[
i
].
Pos
[
k
];
StickyDataIn
[
nexport
].
Vel
[
k
]
=
SphP
[
i
].
VelPred
[
k
];
}
StickyDataIn
[
nexport
].
Mass
=
P
[
i
].
Mass
;
StickyDataIn
[
nexport
].
Hsml
=
SphP
[
i
].
Hsml
;
StickyDataIn
[
nexport
].
StickyMaxID
=
SphP
[
i
].
StickyMaxID
;
StickyDataIn
[
nexport
].
StickyMaxFs
=
SphP
[
i
].
StickyMaxFs
;
StickyDataIn
[
nexport
].
StickyNgb
=
SphP
[
i
].
StickyNgb
;
StickyDataIn
[
nexport
].
ID
=
P
[
i
].
ID
;
StickyDataIn
[
nexport
].
Index
=
i
;
StickyDataIn
[
nexport
].
Task
=
j
;
nexport
++
;
nsend_local
[
j
]
++
;
}
}
}
}
tend
=
second
();
timecomp
+=
timediff
(
tstart
,
tend
);
qsort
(
StickyDataIn
,
nexport
,
sizeof
(
struct
stickydata_in
),
sticky_compare_key
);
for
(
j
=
1
,
noffset
[
0
]
=
0
;
j
<
NTask
;
j
++
)
noffset
[
j
]
=
noffset
[
j
-
1
]
+
nsend_local
[
j
-
1
];
tstart
=
second
();
MPI_Allgather
(
nsend_local
,
NTask
,
MPI_INT
,
nsend
,
NTask
,
MPI_INT
,
MPI_COMM_WORLD
);
tend
=
second
();
timeimbalance
+=
timediff
(
tstart
,
tend
);
/* now do the particles that need to be exported */
for
(
level
=
1
;
level
<
(
1
<<
PTask
);
level
++
)
{
tstart
=
second
();
for
(
j
=
0
;
j
<
NTask
;
j
++
)
nbuffer
[
j
]
=
0
;
for
(
ngrp
=
level
;
ngrp
<
(
1
<<
PTask
);
ngrp
++
)
{
maxfill
=
0
;
for
(
j
=
0
;
j
<
NTask
;
j
++
)
{
if
((
j
^
ngrp
)
<
NTask
)
if
(
maxfill
<
nbuffer
[
j
]
+
nsend
[(
j
^
ngrp
)
*
NTask
+
j
])
maxfill
=
nbuffer
[
j
]
+
nsend
[(
j
^
ngrp
)
*
NTask
+
j
];
}
if
(
maxfill
>=
All
.
BunchSizeSticky
)
break
;
sendTask
=
ThisTask
;
recvTask
=
ThisTask
^
ngrp
;
if
(
recvTask
<
NTask
)
{
if
(
nsend
[
ThisTask
*
NTask
+
recvTask
]
>
0
||
nsend
[
recvTask
*
NTask
+
ThisTask
]
>
0
)
{
/* get the particles */
MPI_Sendrecv
(
&
StickyDataIn
[
noffset
[
recvTask
]],
nsend_local
[
recvTask
]
*
sizeof
(
struct
stickydata_in
),
MPI_BYTE
,
recvTask
,
TAG_HYDRO_A
,
&
StickyDataGet
[
nbuffer
[
ThisTask
]],
nsend
[
recvTask
*
NTask
+
ThisTask
]
*
sizeof
(
struct
stickydata_in
),
MPI_BYTE
,
recvTask
,
TAG_HYDRO_A
,
MPI_COMM_WORLD
,
&
status
);
}
}
for
(
j
=
0
;
j
<
NTask
;
j
++
)
if
((
j
^
ngrp
)
<
NTask
)
nbuffer
[
j
]
+=
nsend
[(
j
^
ngrp
)
*
NTask
+
j
];
}
tend
=
second
();
timecommsumm
+=
timediff
(
tstart
,
tend
);
/* now do the imported particles */
tstart
=
second
();
for
(
j
=
0
;
j
<
nbuffer
[
ThisTask
];
j
++
)
sticky_evaluate
(
j
,
1
,
loop
);
tend
=
second
();
timecomp
+=
timediff
(
tstart
,
tend
);
/* do a block to measure imbalance */
tstart
=
second
();
MPI_Barrier
(
MPI_COMM_WORLD
);
tend
=
second
();
timeimbalance
+=
timediff
(
tstart
,
tend
);
/* get the result */
tstart
=
second
();
for
(
j
=
0
;
j
<
NTask
;
j
++
)
nbuffer
[
j
]
=
0
;
for
(
ngrp
=
level
;
ngrp
<
(
1
<<
PTask
);
ngrp
++
)
{
maxfill
=
0
;
for
(
j
=
0
;
j
<
NTask
;
j
++
)
{
if
((
j
^
ngrp
)
<
NTask
)
if
(
maxfill
<
nbuffer
[
j
]
+
nsend
[(
j
^
ngrp
)
*
NTask
+
j
])
maxfill
=
nbuffer
[
j
]
+
nsend
[(
j
^
ngrp
)
*
NTask
+
j
];
}
if
(
maxfill
>=
All
.
BunchSizeSticky
)
break
;
sendTask
=
ThisTask
;
recvTask
=
ThisTask
^
ngrp
;
if
(
recvTask
<
NTask
)
{
if
(
nsend
[
ThisTask
*
NTask
+
recvTask
]
>
0
||
nsend
[
recvTask
*
NTask
+
ThisTask
]
>
0
)
{
/* send the results */
MPI_Sendrecv
(
&
StickyDataResult
[
nbuffer
[
ThisTask
]],
nsend
[
recvTask
*
NTask
+
ThisTask
]
*
sizeof
(
struct
stickydata_out
),
MPI_BYTE
,
recvTask
,
TAG_HYDRO_B
,
&
StickyDataPartialResult
[
noffset
[
recvTask
]],
nsend_local
[
recvTask
]
*
sizeof
(
struct
stickydata_out
),
MPI_BYTE
,
recvTask
,
TAG_HYDRO_B
,
MPI_COMM_WORLD
,
&
status
);
/* add the result to the particles */
for
(
j
=
0
;
j
<
nsend_local
[
recvTask
];
j
++
)
{
source
=
j
+
noffset
[
recvTask
];
place
=
StickyDataIn
[
source
].
Index
;
if
(
StickyDataPartialResult
[
source
].
StickyMaxFs
>
SphP
[
place
].
StickyMaxFs
)
{
SphP
[
place
].
StickyMaxID
=
StickyDataPartialResult
[
source
].
StickyMaxID
;
SphP
[
place
].
StickyMaxFs
=
StickyDataPartialResult
[
source
].
StickyMaxFs
;
}
if
(
StickyDataPartialResult
[
source
].
StickyNgb
!=-
1
)
/* copy only if met a neighbor */
{
SphP
[
place
].
StickyNgb
=
StickyDataPartialResult
[
source
].
StickyNgb
;
for
(
k
=
0
;
k
<
3
;
k
++
)
SphP
[
place
].
StickyNewVel
[
k
]
=
StickyDataPartialResult
[
source
].
StickyNewVel
[
k
];
}
}
}
}
for
(
j
=
0
;
j
<
NTask
;
j
++
)
if
((
j
^
ngrp
)
<
NTask
)
nbuffer
[
j
]
+=
nsend
[(
j
^
ngrp
)
*
NTask
+
j
];
}
tend
=
second
();
timecommsumm
+=
timediff
(
tstart
,
tend
);
level
=
ngrp
-
1
;
}
MPI_Allgather
(
&
ndone
,
1
,
MPI_INT
,
ndonelist
,
1
,
MPI_INT
,
MPI_COMM_WORLD
);
for
(
j
=
0
;
j
<
NTask
;
j
++
)
ntotleft
-=
ndonelist
[
j
];
}
free
(
ndonelist
);
free
(
nsend
);
free
(
nsend_local
);
free
(
nbuffer
);
free
(
noffset
);
if
(
loop
==
2
)
{
/* do final operations on results */
double
norm_dv
;
double
max_dv
;
int
amax_dv
;
max_dv
=
0
;
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
#ifdef SFR
//if((P[i].Ti_endstep == All.Ti_Current) && (P[i].Type == 0))
#else
//if(P[i].Ti_endstep == All.Ti_Current)
#endif
if
((
SphP
[
i
].
Phase
==
GAS_STICKY
)
&&
(
SphP
[
i
].
StickyFlag
))
{
if
(
loop
==
1
)
{
if
(
SphP
[
i
].
StickyMaxID
!=-
1
)
{
NumColLocal
+=
1
;
}
}
if
(
SphP
[
i
].
StickyNgb
!=-
1
)
{
//printf("(%d) COLLISION : %d -> %d\n",ThisTask,P[i].ID,SphP[i].StickyNgb);
//printf(" : OldV %g NewV %g\n",SphP[i].VelPred[0],SphP[i].StickyNewVel[0]);
//printf(" : OldV %g NewV %g\n",SphP[i].VelPred[1],SphP[i].StickyNewVel[1]);
//printf(" : OldV %g NewV %g\n",SphP[i].VelPred[2],SphP[i].StickyNewVel[2]);
/* new velocities */
for
(
k
=
0
;
k
<
3
;
k
++
)
{
dv
[
k
]
=
SphP
[
i
].
StickyNewVel
[
k
]
-
SphP
[
i
].
VelPred
[
k
];
SphP
[
i
].
VelPred
[
k
]
=
SphP
[
i
].
StickyNewVel
[
k
];
P
[
i
].
Vel
[
k
]
+=
dv
[
k
];
}
norm_dv
=
sqrt
(
dv
[
0
]
*
dv
[
0
]
+
dv
[
1
]
*
dv
[
1
]
+
dv
[
2
]
*
dv
[
2
]);
if
(
norm_dv
>
max_dv
)
{
max_dv
=
norm_dv
;
amax_dv
=
P
[
i
].
ID
;
}
/* set particles as non sticky-active */
SphP
[
i
].
StickyFlag
=
0
;
SphP
[
i
].
StickyTime
=
All
.
Time
+
All
.
StickyIdleTime
;
/* record collision */
NumColLocal
+=
1
;
}
}
printf
(
"(%d) max_dv=%g amax_dv=%d
\n
"
,
ThisTask
,
max_dv
,
amax_dv
);
}
tend
=
second
();
timecomp
+=
timediff
(
tstart
,
tend
);
/* write some statistics */
if
(
loop
==
2
)
{
MPI_Allreduce
(
&
NumColLocal
,
&
NumCol
,
1
,
MPI_INT
,
MPI_SUM
,
MPI_COMM_WORLD
);
MPI_Allreduce
(
&
NumColPotLocal
,
&
NumColPot
,
1
,
MPI_INT
,
MPI_SUM
,
MPI_COMM_WORLD
);
MPI_Allreduce
(
&
NumNoColLocal
,
&
NumNoCol
,
1
,
MPI_INT
,
MPI_SUM
,
MPI_COMM_WORLD
);
if
(
ThisTask
==
0
)
{
fprintf
(
FdSticky
,
"Step %d, Time: %g NumColPot: %d NumCol: %d NumNoCol: %d
\n
"
,
All
.
NumCurrentTiStep
,
All
.
Time
,(
int
)
NumColPot
,(
int
)
NumCol
,(
int
)
NumNoCol
);
fflush
(
FdSticky
);
}
}
if
(
ThisTask
==
0
)
printf
(
"sticky collisions done.
\n
"
);
}
/*!
*
* first sticky
*
*/
void
sticky_evaluate
(
int
target
,
int
mode
,
int
loop
)
{
int
j
,
k
,
n
,
startnode
,
numngb
;
FLOAT
*
pos
,
*
vel
;
FLOAT
mass
;
FLOAT
h_i
;
double
dx
,
dy
,
dz
;
double
dvx
,
dvy
,
dvz
;
double
vdotr
,
vdotr2
,
dv
,
dv2
;
double
vcm
[
3
],
newvel
[
3
];
double
h_i2
;
double
h_j
,
r
,
r2
;
int
phase
=
0
;
int
StickyMaxID
;
int
id
;
double
M1M
,
M2M
,
M12
;
double
dbeta_dv
,
dbeta_dv_er_x
,
dbeta_dv_er_y
,
dbeta_dv_er_z
;
double
dv_beta_t_x
,
dv_beta_t_y
,
dv_beta_t_z
;
double
beta_r
,
beta_t
;
float
fs
,
maxfs
;
int
id_max
;
int
id_ngb
;
beta_r
=
All
.
StickyBetaR
;
beta_t
=
All
.
StickyBetaT
;
if
(
mode
==
0
)
{
pos
=
P
[
target
].
Pos
;
vel
=
SphP
[
target
].
VelPred
;
mass
=
P
[
target
].
Mass
;
h_i
=
SphP
[
target
].
Hsml
;
id
=
P
[
target
].
ID
;
StickyMaxID
=
SphP
[
target
].
StickyMaxID
;
maxfs
=
SphP
[
target
].
StickyMaxFs
;
id_max
=
SphP
[
target
].
StickyMaxID
;
id_ngb
=
SphP
[
target
].
StickyNgb
;
}
else
{
pos
=
StickyDataGet
[
target
].
Pos
;
vel
=
StickyDataGet
[
target
].
Vel
;
mass
=
StickyDataGet
[
target
].
Mass
;
h_i
=
StickyDataGet
[
target
].
Hsml
;
id
=
StickyDataGet
[
target
].
ID
;
StickyMaxID
=
StickyDataGet
[
target
].
StickyMaxID
;
maxfs
=
StickyDataGet
[
target
].
StickyMaxFs
;
id_max
=
StickyDataGet
[
target
].
StickyMaxID
;
id_ngb
=
StickyDataGet
[
target
].
StickyNgb
;
}
h_i2
=
h_i
*
h_i
;
for
(
k
=
0
;
k
<
3
;
k
++
)
newvel
[
k
]
=
vel
[
k
];
/* Now start sticky 1 for this particle */
startnode
=
All
.
MaxPart
;
do
{
numngb
=
ngb_treefind_pairs
(
&
pos
[
0
],
h_i
,
phase
,
&
startnode
);
for
(
n
=
0
;
n
<
numngb
;
n
++
)
{
j
=
Ngblist
[
n
];
dx
=
pos
[
0
]
-
P
[
j
].
Pos
[
0
];
dy
=
pos
[
1
]
-
P
[
j
].
Pos
[
1
];
dz
=
pos
[
2
]
-
P
[
j
].
Pos
[
2
];
#ifdef PERIODIC
/* find the closest image in the given box size */
if
(
dx
>
boxHalf_X
)
dx
-=
boxSize_X
;
if
(
dx
<
-
boxHalf_X
)
dx
+=
boxSize_X
;
if
(
dy
>
boxHalf_Y
)
dy
-=
boxSize_Y
;
if
(
dy
<
-
boxHalf_Y
)
dy
+=
boxSize_Y
;
if
(
dz
>
boxHalf_Z
)
dz
-=
boxSize_Z
;
if
(
dz
<
-
boxHalf_Z
)
dz
+=
boxSize_Z
;
#endif
r2
=
dx
*
dx
+
dy
*
dy
+
dz
*
dz
;
h_j
=
SphP
[
j
].
Hsml
;
if
(
r2
<
h_i2
||
r2
<
h_j
*
h_j
)
{
r
=
sqrt
(
r2
);
if
(
r
>
0
)
{
dvx
=
vel
[
0
]
-
SphP
[
j
].
VelPred
[
0
];
dvy
=
vel
[
1
]
-
SphP
[
j
].
VelPred
[
1
];
dvz
=
vel
[
2
]
-
SphP
[
j
].
VelPred
[
2
];
vdotr
=
dx
*
dvx
+
dy
*
dvy
+
dz
*
dvz
;
vdotr2
=
vdotr
;
dv
=
vdotr2
/
r
;
if
(
loop
==
1
)
{
if
(
vdotr2
<
0
)
/* particles approaches !!! */
{
/* compute new velocities */
M12
=
mass
+
P
[
j
].
Mass
;
M1M
=
mass
/
M12
;
M2M
=
P
[
j
].
Mass
/
M12
;
/* compute velocity of the mass center */
for
(
k
=
0
;
k
<
3
;
k
++
)
vcm
[
k
]
=
(
mass
*
vel
[
k
]
+
P
[
j
].
Mass
*
SphP
[
j
].
VelPred
[
k
])
/
M12
;
dbeta_dv
=
(
beta_r
-
beta_t
)
*
dv
;
dbeta_dv_er_x
=
-
dbeta_dv
*
dx
/
r
;
dbeta_dv_er_y
=
-
dbeta_dv
*
dy
/
r
;
dbeta_dv_er_z
=
-
dbeta_dv
*
dz
/
r
;
dv_beta_t_x
=
dvx
*
beta_t
;
dv_beta_t_y
=
dvy
*
beta_t
;
dv_beta_t_z
=
dvz
*
beta_t
;
newvel
[
0
]
=
M2M
*
(
+
dv_beta_t_x
-
dbeta_dv_er_x
)
+
vcm
[
0
];
newvel
[
1
]
=
M2M
*
(
+
dv_beta_t_y
-
dbeta_dv_er_y
)
+
vcm
[
1
];
newvel
[
2
]
=
M2M
*
(
+
dv_beta_t_z
-
dbeta_dv_er_z
)
+
vcm
[
2
];
/* check that dv is not too large */
dv2
=
(
newvel
[
0
]
-
vel
[
0
])
*
(
newvel
[
0
]
-
vel
[
0
])
+
(
newvel
[
1
]
-
vel
[
1
])
*
(
newvel
[
1
]
-
vel
[
1
])
+
(
newvel
[
2
]
-
vel
[
2
])
*
(
newvel
[
2
]
-
vel
[
2
]);
if
(
dv2
<
All
.
StickyMaxVelocity
*
All
.
StickyMaxVelocity
)
{
fs
=
r
/
(
(
0.5
*
(
h_i
+
h_j
)
+
r
)
*
(
0.5
*
(
h_i
+
h_j
)
+
r
));
if
(
fs
>
maxfs
)
{
id_max
=
P
[
j
].
ID
;
maxfs
=
fs
;
}
}
}
}
else
{
if
((
SphP
[
j
].
StickyMaxID
==
id
)
&&
(
StickyMaxID
==
P
[
j
].
ID
))
{
id_ngb
=
P
[
j
].
ID
;
/* compute new velocities */
M12
=
mass
+
P
[
j
].
Mass
;
M1M
=
mass
/
M12
;
M2M
=
P
[
j
].
Mass
/
M12
;
/* compute velocity of the mass center */
for
(
k
=
0
;
k
<
3
;
k
++
)
vcm
[
k
]
=
(
mass
*
vel
[
k
]
+
P
[
j
].
Mass
*
SphP
[
j
].
VelPred
[
k
])
/
M12
;
dbeta_dv
=
(
beta_r
-
beta_t
)
*
dv
;
dbeta_dv_er_x
=
-
dbeta_dv
*
dx
/
r
;
dbeta_dv_er_y
=
-
dbeta_dv
*
dy
/
r
;
dbeta_dv_er_z
=
-
dbeta_dv
*
dz
/
r
;
dv_beta_t_x
=
dvx
*
beta_t
;
dv_beta_t_y
=
dvy
*
beta_t
;
dv_beta_t_z
=
dvz
*
beta_t
;
newvel
[
0
]
=
M2M
*
(
+
dv_beta_t_x
-
dbeta_dv_er_x
)
+
vcm
[
0
];
newvel
[
1
]
=
M2M
*
(
+
dv_beta_t_y
-
dbeta_dv_er_y
)
+
vcm
[
1
];
newvel
[
2
]
=
M2M
*
(
+
dv_beta_t_z
-
dbeta_dv_er_z
)
+
vcm
[
2
];
}
}
}
}
}
}
while
(
startnode
>=
0
);
/* Now collect the result at the right place */
if
(
mode
==
0
)
{
if
(
loop
==
1
)
{
SphP
[
target
].
StickyMaxID
=
id_max
;
SphP
[
target
].
StickyMaxFs
=
maxfs
;
}
SphP
[
target
].
StickyNgb
=
id_ngb
;
for
(
k
=
0
;
k
<
3
;
k
++
)
SphP
[
target
].
StickyNewVel
[
k
]
=
newvel
[
k
];
}
else
{
if
(
loop
==
1
)
{
StickyDataResult
[
target
].
StickyMaxID
=
id_max
;
StickyDataResult
[
target
].
StickyMaxFs
=
maxfs
;
}
StickyDataResult
[
target
].
StickyNgb
=
id_ngb
;
for
(
k
=
0
;
k
<
3
;
k
++
)
StickyDataResult
[
target
].
StickyNewVel
[
k
]
=
newvel
[
k
];
}
}
/*! This is a comparison kernel for a sort routine, which is used to group
* particles that are going to be exported to the same CPU.
*/
int
sticky_compare_key
(
const
void
*
a
,
const
void
*
b
)
{
if
(((
struct
stickydata_in
*
)
a
)
->
Task
<
(((
struct
stickydata_in
*
)
b
)
->
Task
))
return
-
1
;
if
(((
struct
stickydata_in
*
)
a
)
->
Task
>
(((
struct
stickydata_in
*
)
b
)
->
Task
))
return
+
1
;
return
0
;
}
#endif
Event Timeline
Log In to Comment