Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77958902
stars_density.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
Sat, Aug 17, 11:52
Size
19 KB
Mime Type
text/x-c
Expires
Mon, Aug 19, 11:52 (2 d)
Engine
blob
Format
Raw Data
Handle
19956550
Attached To
rGEAR Gear
stars_density.c
View Options
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
/*! \file density.c
* \brief SPH density computation and smoothing length determination
*
* This file contains the "first SPH loop", where the SPH densities and
* some auxiliary quantities are computed. If the number of neighbours
* obtained falls outside the target range, the correct smoothing
* length is determined iteratively, if needed.
*/
#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 CHIMIE
/*! This function computes the local density for each active SPH particle,
* the number of neighbours in the current smoothing radius, and the
* divergence and curl of the velocity field. The pressure is updated as
* well. If a particle with its smoothing region is fully inside the
* local domain, it is not exported to the other processors. The function
* also detects particles that have a number of neighbours outside the
* allowed tolerance range. For these particles, the smoothing length is
* adjusted accordingly, and the density computation is executed again.
* Note that the smoothing length is not allowed to fall below the lower
* bound set by MinGasHsml.
*/
void
stars_density
(
void
)
{
long
long
ntot
,
ntotleft
;
int
*
noffset
,
*
nbuffer
,
*
nsend
,
*
nsend_local
,
*
numlist
,
*
ndonelist
;
int
i
,
j
,
n
,
m
,
ndone
,
npleft
,
maxfill
,
source
,
iter
=
0
;
int
level
,
ngrp
,
sendTask
,
recvTask
,
place
,
nexport
;
double
tstart
,
tend
,
tstart_ngb
=
0
,
tend_ngb
=
0
;
double
sumt
,
sumcomm
,
timengb
,
sumtimengb
;
double
timecomp
=
0
,
timeimbalance
=
0
,
timecommsumm
=
0
,
sumimbalance
;
MPI_Status
status
;
#ifdef DETAILED_CPU_OUTPUT_IN_STARS_DENSITY
double
*
timengblist
;
double
*
timecomplist
;
double
*
timecommsummlist
;
double
*
timeimbalancelist
;
#endif
#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
if
(
ThisTask
==
0
)
printf
(
"Start star density computation.
\n
"
);
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
);
/* `NumStUpdate' gives the number of particles on this processor that want a chimie computation */
for
(
n
=
N_gas
,
NumStUpdate
=
0
;
n
<
N_gas
+
N_stars
;
n
++
)
{
m
=
P
[
n
].
StPIdx
;
StP
[
m
].
Left
=
StP
[
m
].
Right
=
0
;
#ifdef AVOIDNUMNGBPROBLEM
StP
[
m
].
OldNumNgb
=
-
1
;
#endif
if
(
P
[
n
].
Ti_endstep
==
All
.
Ti_Current
)
NumStUpdate
++
;
}
numlist
=
malloc
(
NTask
*
sizeof
(
int
)
*
NTask
);
MPI_Allgather
(
&
NumStUpdate
,
1
,
MPI_INT
,
numlist
,
1
,
MPI_INT
,
MPI_COMM_WORLD
);
for
(
i
=
0
,
ntot
=
0
;
i
<
NTask
;
i
++
)
ntot
+=
numlist
[
i
];
free
(
numlist
);
/* we will repeat the whole thing for those particles where we didn't
* find enough neighbours
*/
do
{
i
=
N_gas
;
/* beginn with the first star */
ntotleft
=
ntot
;
/* particles left for all tasks together */
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
+
N_stars
&&
nexport
<
All
.
BunchSizeDensity
-
NTask
;
i
++
)
if
(
P
[
i
].
Ti_endstep
==
All
.
Ti_Current
)
{
if
(
P
[
i
].
Type
!=
ST
)
{
printf
(
"P[i].Type != ST, we better stop.
\n
"
);
printf
(
"N_gas=%d (type=%d) i=%d (type=%d)
\n
"
,
N_gas
,
P
[
N_gas
].
Type
,
i
,
P
[
i
].
Type
);
printf
(
"Please, check that you do not use PEANOHILBERT
\n
"
);
endrun
(
666001
);
}
ndone
++
;
m
=
P
[
i
].
StPIdx
;
for
(
j
=
0
;
j
<
NTask
;
j
++
)
Exportflag
[
j
]
=
0
;
stars_density_evaluate
(
i
,
0
);
for
(
j
=
0
;
j
<
NTask
;
j
++
)
{
if
(
Exportflag
[
j
])
{
StarsDensDataIn
[
nexport
].
Pos
[
0
]
=
P
[
i
].
Pos
[
0
];
StarsDensDataIn
[
nexport
].
Pos
[
1
]
=
P
[
i
].
Pos
[
1
];
StarsDensDataIn
[
nexport
].
Pos
[
2
]
=
P
[
i
].
Pos
[
2
];
StarsDensDataIn
[
nexport
].
Hsml
=
StP
[
m
].
Hsml
;
StarsDensDataIn
[
nexport
].
Index
=
i
;
StarsDensDataIn
[
nexport
].
Task
=
j
;
nexport
++
;
nsend_local
[
j
]
++
;
}
}
}
tend
=
second
();
timecomp
+=
timediff
(
tstart
,
tend
);
qsort
(
StarsDensDataIn
,
nexport
,
sizeof
(
struct
starsdensdata_in
),
stars_dens_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
.
BunchSizeDensity
)
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
(
&
StarsDensDataIn
[
noffset
[
recvTask
]],
nsend_local
[
recvTask
]
*
sizeof
(
struct
starsdensdata_in
),
MPI_BYTE
,
recvTask
,
TAG_DENS_A
,
&
StarsDensDataGet
[
nbuffer
[
ThisTask
]],
nsend
[
recvTask
*
NTask
+
ThisTask
]
*
sizeof
(
struct
starsdensdata_in
),
MPI_BYTE
,
recvTask
,
TAG_DENS_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
);
tstart
=
second
();
for
(
j
=
0
;
j
<
nbuffer
[
ThisTask
];
j
++
)
stars_density_evaluate
(
j
,
1
);
tend
=
second
();
timecomp
+=
timediff
(
tstart
,
tend
);
/* do a block to explicitly 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
.
BunchSizeDensity
)
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
(
&
StarsDensDataResult
[
nbuffer
[
ThisTask
]],
nsend
[
recvTask
*
NTask
+
ThisTask
]
*
sizeof
(
struct
starsdensdata_out
),
MPI_BYTE
,
recvTask
,
TAG_DENS_B
,
&
StarsDensDataPartialResult
[
noffset
[
recvTask
]],
nsend_local
[
recvTask
]
*
sizeof
(
struct
starsdensdata_out
),
MPI_BYTE
,
recvTask
,
TAG_DENS_B
,
MPI_COMM_WORLD
,
&
status
);
/* add the result to the particles */
for
(
j
=
0
;
j
<
nsend_local
[
recvTask
];
j
++
)
{
source
=
j
+
noffset
[
recvTask
];
place
=
P
[
StarsDensDataIn
[
source
].
Index
].
StPIdx
;
StP
[
place
].
NumNgb
+=
StarsDensDataPartialResult
[
source
].
Ngb
;
StP
[
place
].
Density
+=
StarsDensDataPartialResult
[
source
].
Rho
;
StP
[
place
].
Volume
+=
StarsDensDataPartialResult
[
source
].
Volume
;
#ifdef CHIMIE_KINETIC_FEEDBACK
StP
[
place
].
NgbMass
+=
StarsDensDataPartialResult
[
source
].
NgbMass
;
#endif
StP
[
place
].
DhsmlDensityFactor
+=
StarsDensDataPartialResult
[
source
].
DhsmlDensity
;
}
}
}
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
];
}
/* do final operations on results */
tstart
=
second
();
for
(
i
=
N_gas
,
npleft
=
0
;
i
<
N_gas
+
N_stars
;
i
++
)
{
if
((
P
[
i
].
Ti_endstep
==
All
.
Ti_Current
))
{
m
=
P
[
i
].
StPIdx
;
StP
[
m
].
DhsmlDensityFactor
=
1
/
(
1
+
StP
[
m
].
Hsml
*
StP
[
m
].
DhsmlDensityFactor
/
(
NUMDIMS
*
StP
[
m
].
Density
));
/* now check whether we had enough neighbours */
if
(
StP
[
m
].
NumNgb
<
(
All
.
DesNumNgb
-
All
.
MaxNumNgbDeviation
)
||
(
StP
[
m
].
NumNgb
>
(
All
.
DesNumNgb
+
All
.
MaxNumNgbDeviation
)
&&
StP
[
m
].
Hsml
>
(
1.01
*
All
.
MinGasHsml
)))
{
#ifdef AVOIDNUMNGBPROBLEM
// if((StP[m].NumNgb>StP[m].OldNumNgb-All.MaxNumNgbDeviation/10000.)
// &&(StP[m].NumNgb<StP[m].OldNumNgb+All.MaxNumNgbDeviation/10000.))
if
(
StP
[
m
].
NumNgb
==
StP
[
m
].
OldNumNgb
)
{
P
[
i
].
Ti_endstep
=
-
P
[
i
].
Ti_endstep
-
1
;
printf
(
"ID=%d NumNgb=%g OldNumNgb=%g
\n
"
,
P
[
i
].
ID
,
StP
[
m
].
NumNgb
,
StP
[
m
].
OldNumNgb
);
continue
;
}
StP
[
m
].
OldNumNgb
=
StP
[
m
].
NumNgb
;
#endif
/* need to redo this particle */
npleft
++
;
if
(
StP
[
m
].
Left
>
0
&&
StP
[
m
].
Right
>
0
)
if
((
StP
[
m
].
Right
-
StP
[
m
].
Left
)
<
1.0e-3
*
StP
[
m
].
Left
)
{
/* this one should be ok */
npleft
--
;
P
[
i
].
Ti_endstep
=
-
P
[
i
].
Ti_endstep
-
1
;
/* Mark as inactive */
continue
;
}
if
(
StP
[
m
].
NumNgb
<
(
All
.
DesNumNgb
-
All
.
MaxNumNgbDeviation
))
StP
[
m
].
Left
=
dmax
(
StP
[
m
].
Hsml
,
StP
[
m
].
Left
);
else
{
if
(
StP
[
m
].
Right
!=
0
)
{
if
(
StP
[
m
].
Hsml
<
StP
[
m
].
Right
)
StP
[
m
].
Right
=
StP
[
m
].
Hsml
;
}
else
StP
[
m
].
Right
=
StP
[
m
].
Hsml
;
}
if
(
iter
>=
MAXITER
-
10
)
{
printf
(
"i=%d task=%d ID=%d Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g
\n
pos=(%g|%g|%g)
\n
"
,
i
,
ThisTask
,
(
int
)
P
[
i
].
ID
,
StP
[
m
].
Hsml
,
StP
[
m
].
Left
,
StP
[
m
].
Right
,
(
float
)
StP
[
m
].
NumNgb
,
StP
[
m
].
Right
-
StP
[
m
].
Left
,
P
[
i
].
Pos
[
0
],
P
[
i
].
Pos
[
1
],
P
[
i
].
Pos
[
2
]);
fflush
(
stdout
);
}
if
(
StP
[
m
].
Right
>
0
&&
StP
[
m
].
Left
>
0
)
StP
[
m
].
Hsml
=
pow
(
0.5
*
(
pow
(
StP
[
m
].
Left
,
3
)
+
pow
(
StP
[
m
].
Right
,
3
)),
1.0
/
3
);
else
{
if
(
StP
[
m
].
Right
==
0
&&
StP
[
m
].
Left
==
0
)
{
printf
(
"i=%d task=%d ID=%d Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g
\n
pos=(%g|%g|%g)
\n
"
,
i
,
ThisTask
,
(
int
)
P
[
i
].
ID
,
StP
[
m
].
Hsml
,
StP
[
m
].
Left
,
StP
[
m
].
Right
,
(
float
)
StP
[
m
].
NumNgb
,
StP
[
m
].
Right
-
StP
[
m
].
Left
,
P
[
i
].
Pos
[
0
],
P
[
i
].
Pos
[
1
],
P
[
i
].
Pos
[
2
]);
fflush
(
stdout
);
endrun
(
8188
);
/* can't occur */
}
if
(
StP
[
m
].
Right
==
0
&&
StP
[
m
].
Left
>
0
)
{
if
(
P
[
i
].
Type
==
ST
&&
fabs
(
StP
[
m
].
NumNgb
-
All
.
DesNumNgb
)
<
0.5
*
All
.
DesNumNgb
)
{
StP
[
m
].
Hsml
*=
1
-
(
StP
[
m
].
NumNgb
-
All
.
DesNumNgb
)
/
(
NUMDIMS
*
StP
[
m
].
NumNgb
)
*
StP
[
m
].
DhsmlDensityFactor
;
}
else
StP
[
m
].
Hsml
*=
1.26
;
}
if
(
StP
[
m
].
Right
>
0
&&
StP
[
m
].
Left
==
0
)
{
if
(
P
[
i
].
Type
==
ST
&&
fabs
(
StP
[
m
].
NumNgb
-
All
.
DesNumNgb
)
<
0.5
*
All
.
DesNumNgb
)
{
StP
[
m
].
Hsml
*=
1
-
(
StP
[
m
].
NumNgb
-
All
.
DesNumNgb
)
/
(
NUMDIMS
*
StP
[
m
].
NumNgb
)
*
StP
[
m
].
DhsmlDensityFactor
;
}
else
StP
[
m
].
Hsml
/=
1.26
;
}
}
if
(
StP
[
m
].
Hsml
<
All
.
MinGasHsml
)
StP
[
m
].
Hsml
=
All
.
MinGasHsml
;
}
else
P
[
i
].
Ti_endstep
=
-
P
[
i
].
Ti_endstep
-
1
;
/* Mark as inactive */
}
}
tend
=
second
();
timecomp
+=
timediff
(
tstart
,
tend
);
numlist
=
malloc
(
NTask
*
sizeof
(
int
)
*
NTask
);
MPI_Allgather
(
&
npleft
,
1
,
MPI_INT
,
numlist
,
1
,
MPI_INT
,
MPI_COMM_WORLD
);
for
(
i
=
0
,
ntot
=
0
;
i
<
NTask
;
i
++
)
ntot
+=
numlist
[
i
];
free
(
numlist
);
if
(
ntot
>
0
)
{
if
(
iter
==
0
)
tstart_ngb
=
second
();
iter
++
;
if
(
iter
>
0
&&
ThisTask
==
0
)
{
printf
(
"ngb iteration %d: need to repeat for %d%09d particles.
\n
"
,
iter
,
(
int
)
(
ntot
/
1000000000
),
(
int
)
(
ntot
%
1000000000
));
fflush
(
stdout
);
}
if
(
iter
>
MAXITER
)
{
printf
(
"failed to converge in neighbour iteration in density()
\n
"
);
fflush
(
stdout
);
endrun
(
1155
);
}
}
else
tend_ngb
=
second
();
}
while
(
ntot
>
0
);
/* mark as active again */
for
(
i
=
0
;
i
<
NumPart
;
i
++
)
if
(
P
[
i
].
Ti_endstep
<
0
)
P
[
i
].
Ti_endstep
=
-
P
[
i
].
Ti_endstep
-
1
;
free
(
ndonelist
);
free
(
nsend
);
free
(
nsend_local
);
free
(
nbuffer
);
free
(
noffset
);
/* collect some timing information */
if
(
iter
>
0
)
timengb
=
timediff
(
tstart_ngb
,
tend_ngb
);
else
timengb
=
0
;
MPI_Reduce
(
&
timengb
,
&
sumtimengb
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
0
,
MPI_COMM_WORLD
);
MPI_Reduce
(
&
timecomp
,
&
sumt
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
0
,
MPI_COMM_WORLD
);
MPI_Reduce
(
&
timecommsumm
,
&
sumcomm
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
0
,
MPI_COMM_WORLD
);
MPI_Reduce
(
&
timeimbalance
,
&
sumimbalance
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
0
,
MPI_COMM_WORLD
);
#ifdef DETAILED_CPU_OUTPUT_IN_STARS_DENSITY
numlist
=
malloc
(
sizeof
(
int
)
*
NTask
);
timengblist
=
malloc
(
sizeof
(
double
)
*
NTask
);
timecomplist
=
malloc
(
sizeof
(
double
)
*
NTask
);
timecommsummlist
=
malloc
(
sizeof
(
double
)
*
NTask
);
timeimbalancelist
=
malloc
(
sizeof
(
double
)
*
NTask
);
MPI_Gather
(
&
NumStUpdate
,
1
,
MPI_INT
,
numlist
,
1
,
MPI_INT
,
0
,
MPI_COMM_WORLD
);
MPI_Gather
(
&
timengb
,
1
,
MPI_DOUBLE
,
timengblist
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
MPI_Gather
(
&
timecomp
,
1
,
MPI_DOUBLE
,
timecomplist
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
MPI_Gather
(
&
timecommsumm
,
1
,
MPI_DOUBLE
,
timecommsummlist
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
MPI_Gather
(
&
timeimbalance
,
1
,
MPI_DOUBLE
,
timeimbalancelist
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
if
(
ThisTask
==
0
)
{
fprintf
(
FdTimings
,
"
\n
stars density
\n\n
"
);
fprintf
(
FdTimings
,
"Nupdate "
);
for
(
i
=
0
;
i
<
NTask
;
i
++
)
fprintf
(
FdTimings
,
"%12d "
,
numlist
[
i
]);
fprintf
(
FdTimings
,
"
\n
"
);
fprintf
(
FdTimings
,
"timengb "
);
for
(
i
=
0
;
i
<
NTask
;
i
++
)
fprintf
(
FdTimings
,
"%12g "
,
timengblist
[
i
]);
fprintf
(
FdTimings
,
"
\n
"
);
fprintf
(
FdTimings
,
"timecomp "
);
for
(
i
=
0
;
i
<
NTask
;
i
++
)
fprintf
(
FdTimings
,
"%12g "
,
timecomplist
[
i
]);
fprintf
(
FdTimings
,
"
\n
"
);
fprintf
(
FdTimings
,
"timecommsumm "
);
for
(
i
=
0
;
i
<
NTask
;
i
++
)
fprintf
(
FdTimings
,
"%12g "
,
timecommsummlist
[
i
]);
fprintf
(
FdTimings
,
"
\n
"
);
fprintf
(
FdTimings
,
"timeimbalance "
);
for
(
i
=
0
;
i
<
NTask
;
i
++
)
fprintf
(
FdTimings
,
"%12g "
,
timeimbalancelist
[
i
]);
fprintf
(
FdTimings
,
"
\n
"
);
fprintf
(
FdTimings
,
"
\n
"
);
}
free
(
timeimbalancelist
);
free
(
timecommsummlist
);
free
(
timecomplist
);
free
(
timengblist
);
free
(
numlist
);
#endif
if
(
ThisTask
==
0
)
{
All
.
CPU_ChimieDensCompWalk
+=
sumt
/
NTask
;
All
.
CPU_ChimieDensCommSumm
+=
sumcomm
/
NTask
;
All
.
CPU_ChimieDensImbalance
+=
sumimbalance
/
NTask
;
All
.
CPU_ChimieDensEnsureNgb
+=
sumtimengb
/
NTask
;
}
if
(
ThisTask
==
0
)
printf
(
"Stars density computation done.
\n
"
);
}
/*! This function represents the core of the SPH density computation. The
* target particle may either be local, or reside in the communication
* buffer.
*/
void
stars_density_evaluate
(
int
target
,
int
mode
)
{
int
j
,
n
,
startnode
,
numngb
,
numngb_inbox
;
double
h
,
h2
,
hinv
,
hinv3
,
hinv4
;
double
rho
,
volume
,
wk
,
dwk
;
double
dx
,
dy
,
dz
,
r
,
r2
,
u
,
mass_j
,
rho_j
;
double
weighted_numngb
,
dhsmlrho
;
int
target_stp
;
FLOAT
*
pos
;
#ifdef CHIMIE_KINETIC_FEEDBACK
double
ngbmass
;
#endif
if
(
mode
==
0
)
{
pos
=
P
[
target
].
Pos
;
target_stp
=
P
[
target
].
StPIdx
;
h
=
StP
[
target_stp
].
Hsml
;
}
else
{
pos
=
StarsDensDataGet
[
target
].
Pos
;
h
=
StarsDensDataGet
[
target
].
Hsml
;
}
h2
=
h
*
h
;
hinv
=
1.0
/
h
;
#ifndef TWODIMS
hinv3
=
hinv
*
hinv
*
hinv
;
#else
hinv3
=
hinv
*
hinv
/
boxSize_Z
;
#endif
hinv4
=
hinv3
*
hinv
;
rho
=
0
;
volume
=
0
;
weighted_numngb
=
0
;
dhsmlrho
=
0
;
#ifdef CHIMIE_KINETIC_FEEDBACK
ngbmass
=
0
;
#endif
startnode
=
All
.
MaxPart
;
numngb
=
0
;
do
{
numngb_inbox
=
ngb_treefind_variable_for_chimie
(
&
pos
[
0
],
h
,
&
startnode
);
for
(
n
=
0
;
n
<
numngb_inbox
;
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
/* now 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
;
if
(
r2
<
h2
)
{
numngb
++
;
r
=
sqrt
(
r2
);
u
=
r
*
hinv
;
if
(
u
<
0.5
)
{
wk
=
hinv3
*
(
KERNEL_COEFF_1
+
KERNEL_COEFF_2
*
(
u
-
1
)
*
u
*
u
);
dwk
=
hinv4
*
u
*
(
KERNEL_COEFF_3
*
u
-
KERNEL_COEFF_4
);
}
else
{
wk
=
hinv3
*
KERNEL_COEFF_5
*
(
1.0
-
u
)
*
(
1.0
-
u
)
*
(
1.0
-
u
);
dwk
=
hinv4
*
KERNEL_COEFF_6
*
(
1.0
-
u
)
*
(
1.0
-
u
);
}
mass_j
=
P
[
j
].
Mass
;
rho_j
=
SphP
[
j
].
Density
;
rho
+=
mass_j
*
wk
;
volume
+=
mass_j
/
rho_j
*
wk
;
weighted_numngb
+=
NORM_COEFF
*
wk
/
hinv3
;
dhsmlrho
+=
-
mass_j
*
(
NUMDIMS
*
hinv
*
wk
+
u
*
dwk
);
#ifdef CHIMIE_KINETIC_FEEDBACK
ngbmass
+=
mass_j
;
#endif
}
}
}
while
(
startnode
>=
0
);
if
(
mode
==
0
)
{
StP
[
target_stp
].
NumNgb
=
weighted_numngb
;
StP
[
target_stp
].
Density
=
rho
;
StP
[
target_stp
].
Volume
=
volume
;
StP
[
target_stp
].
DhsmlDensityFactor
=
dhsmlrho
;
#ifdef CHIMIE_KINETIC_FEEDBACK
StP
[
target_stp
].
NgbMass
=
ngbmass
;
#endif
}
else
{
StarsDensDataResult
[
target
].
Rho
=
rho
;
StarsDensDataResult
[
target
].
Volume
=
volume
;
StarsDensDataResult
[
target
].
Ngb
=
weighted_numngb
;
StarsDensDataResult
[
target
].
DhsmlDensity
=
dhsmlrho
;
#ifdef CHIMIE_KINETIC_FEEDBACK
StarsDensDataResult
[
target
].
NgbMass
=
ngbmass
;
#endif
}
}
/*! This routine is a comparison kernel used in a sort routine to group
* particles that are exported to the same processor.
*/
int
stars_dens_compare_key
(
const
void
*
a
,
const
void
*
b
)
{
if
(((
struct
starsdensdata_in
*
)
a
)
->
Task
<
(((
struct
starsdensdata_in
*
)
b
)
->
Task
))
return
-
1
;
if
(((
struct
starsdensdata_in
*
)
a
)
->
Task
>
(((
struct
starsdensdata_in
*
)
b
)
->
Task
))
return
+
1
;
return
0
;
}
#endif
Event Timeline
Log In to Comment