Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88722568
predict.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
Sun, Oct 20, 08:56
Size
5 KB
Mime Type
text/x-c
Expires
Tue, Oct 22, 08:56 (2 d)
Engine
blob
Format
Raw Data
Handle
21812013
Attached To
rGEAR Gear
predict.c
View Options
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_math.h>
#include "allvars.h"
#include "proto.h"
/*! \file predict.c
* \brief drift particles by a small time interval
*
* This function contains code to implement a drift operation on all the
* particles, which represents one part of the leapfrog integration scheme.
*/
/*! This function drifts all particles from the current time to the future:
* time0 - > time1
*
* If there is no explicit tree construction in the following timestep, the
* tree nodes are also drifted and updated accordingly. Note: For periodic
* boundary conditions, the mapping of coordinates onto the interval
* [0,All.BoxSize] is only done before the domain decomposition, or for
* outputs to snapshot files. This simplifies dynamic tree updates, and
* allows the domain decomposition to be carried out only every once in a
* while.
*/
void
move_particles
(
int
time0
,
int
time1
)
{
int
i
,
j
;
double
dt_drift
,
dt_gravkick
,
dt_hydrokick
,
dt_entr
;
double
t0
,
t1
;
#ifdef METAL_DIFFUSION
double
eAdt
;
#endif
t0
=
second
();
if
(
All
.
ComovingIntegrationOn
)
{
dt_drift
=
get_drift_factor
(
time0
,
time1
);
dt_gravkick
=
get_gravkick_factor
(
time0
,
time1
);
dt_hydrokick
=
get_hydrokick_factor
(
time0
,
time1
);
}
else
{
dt_drift
=
dt_gravkick
=
dt_hydrokick
=
(
time1
-
time0
)
*
All
.
Timebase_interval
;
}
for
(
i
=
0
;
i
<
NumPart
;
i
++
)
{
for
(
j
=
0
;
j
<
3
;
j
++
)
P
[
i
].
Pos
[
j
]
+=
P
[
i
].
Vel
[
j
]
*
dt_drift
;
if
(
P
[
i
].
Type
==
0
)
{
#ifdef PMGRID
for
(
j
=
0
;
j
<
3
;
j
++
)
{
SphP
[
i
].
VelPred
[
j
]
+=
(
P
[
i
].
GravAccel
[
j
]
+
P
[
i
].
GravPM
[
j
])
*
dt_gravkick
+
SphP
[
i
].
HydroAccel
[
j
]
*
dt_hydrokick
;
}
#else
for
(
j
=
0
;
j
<
3
;
j
++
)
{
SphP
[
i
].
VelPred
[
j
]
+=
P
[
i
].
GravAccel
[
j
]
*
dt_gravkick
+
SphP
[
i
].
HydroAccel
[
j
]
*
dt_hydrokick
;
}
#endif
#ifdef DISSIPATION_FORCES
for
(
j
=
0
;
j
<
3
;
j
++
)
{
SphP
[
i
].
VelPred
[
j
]
+=
SphP
[
i
].
DissipationForcesAccel
[
j
]
*
dt_hydrokick
;
}
#endif
SphP
[
i
].
Density
*=
exp
(
-
SphP
[
i
].
DivVel
*
dt_drift
);
SphP
[
i
].
Hsml
*=
exp
(
0.333333333333
*
SphP
[
i
].
DivVel
*
dt_drift
);
if
(
SphP
[
i
].
Hsml
<
All
.
MinGasHsml
)
SphP
[
i
].
Hsml
=
All
.
MinGasHsml
;
dt_entr
=
(
time1
-
(
P
[
i
].
Ti_begstep
+
P
[
i
].
Ti_endstep
)
/
2
)
*
All
.
Timebase_interval
;
#ifdef DENSITY_INDEPENDENT_SPH
SphP
[
i
].
EgyWtDensity
*=
exp
(
-
SphP
[
i
].
DivVel
*
dt_drift
);
SphP
[
i
].
EntVarPred
=
pow
(
SphP
[
i
].
Entropy
+
SphP
[
i
].
DtEntropy
*
dt_entr
,
1
/
GAMMA
);
SphP
[
i
].
Pressure
=
(
SphP
[
i
].
Entropy
+
SphP
[
i
].
DtEntropy
*
dt_entr
)
*
pow
(
SphP
[
i
].
EgyWtDensity
,
GAMMA
);
#else
SphP
[
i
].
Pressure
=
(
SphP
[
i
].
Entropy
+
SphP
[
i
].
DtEntropy
*
dt_entr
)
*
pow
(
SphP
[
i
].
Density
,
GAMMA
);
#endif
#ifdef ENTROPYPRED
SphP
[
i
].
EntropyPred
=
(
SphP
[
i
].
Entropy
+
SphP
[
i
].
DtEntropy
*
dt_entr
);
#endif
#ifdef CHECK_ENTROPY_SIGN
if
((
SphP
[
i
].
EntropyPred
<
0
)
||
(
SphP
[
i
].
Entropy
<
0
))
{
printf
(
"
\n
task=%d: EntropyPred less than zero in move_particles !
\n
"
,
ThisTask
);
printf
(
"ID=%d Entropy=%g EntropyPred=%g DtEntropy=%g dt_entr=%g
\n
"
,
P
[
i
].
ID
,
SphP
[
i
].
Entropy
,
SphP
[
i
].
EntropyPred
,
SphP
[
i
].
DtEntropy
,
dt_entr
);
fflush
(
stdout
);
endrun
(
333021
);
}
#endif
#ifdef NO_NEGATIVE_PRESSURE
if
(
SphP
[
i
].
Pressure
<
0
)
{
printf
(
"
\n
task=%d: pressure less than zero in move_particles !
\n
"
,
ThisTask
);
printf
(
"ID=%d Entropy=%g DtEntropy*dt=%g Density=%g DtEntropy=%g dt=%g
\n
"
,
P
[
i
].
ID
,
SphP
[
i
].
Entropy
,
SphP
[
i
].
DtEntropy
*
dt_entr
,
SphP
[
i
].
Density
,
SphP
[
i
].
DtEntropy
,
dt_entr
);
fflush
(
stdout
);
endrun
(
333022
);
}
#endif
/***********************************************************/
/* compute art visc coeff */
/***********************************************************/
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
move_art_visc
(
i
,
dt_drift
);
#endif
/***********************************************************/
/* compute diffusion */
/***********************************************************/
#ifdef METAL_DIFFUSION
if
(
SphP
[
i
].
MetalDiffusionA
!=
0
)
for
(
j
=
0
;
j
<
NELEMENTS
;
j
++
)
if
(
SphP
[
i
].
MetalDiffusionA
!=
0
)
{
eAdt
=
exp
(
SphP
[
i
].
MetalDiffusionA
*
dt_drift
);
SphP
[
i
].
Metal
[
j
]
=
SphP
[
i
].
Metal
[
j
]
*
eAdt
+
SphP
[
i
].
MetalDiffusionB
[
j
]
/
SphP
[
i
].
MetalDiffusionA
*
(
1.0
-
eAdt
)
;
}
#endif
}
}
/* if domain-decomp and tree are not going to be reconstructed, update dynamically. */
if
(
All
.
NumForcesSinceLastDomainDecomp
<
All
.
TotNumPart
*
All
.
TreeDomainUpdateFrequency
)
{
for
(
i
=
0
;
i
<
Numnodestree
;
i
++
)
for
(
j
=
0
;
j
<
3
;
j
++
)
Nodes
[
All
.
MaxPart
+
i
].
u
.
d
.
s
[
j
]
+=
Extnodes
[
All
.
MaxPart
+
i
].
vs
[
j
]
*
dt_drift
;
force_update_len
();
force_update_pseudoparticles
();
}
t1
=
second
();
All
.
CPU_Predict
+=
timediff
(
t0
,
t1
);
}
/*! This function makes sure that all particle coordinates (Pos) are
* periodically mapped onto the interval [0, BoxSize]. After this function
* has been called, a new domain decomposition should be done, which will
* also force a new tree construction.
*/
#ifdef PERIODIC
void
do_box_wrapping
(
void
)
{
int
i
,
j
;
double
boxsize
[
3
];
for
(
j
=
0
;
j
<
3
;
j
++
)
boxsize
[
j
]
=
All
.
BoxSize
;
#ifdef LONG_X
boxsize
[
0
]
*=
LONG_X
;
#endif
#ifdef LONG_Y
boxsize
[
1
]
*=
LONG_Y
;
#endif
#ifdef LONG_Z
boxsize
[
2
]
*=
LONG_Z
;
#endif
for
(
i
=
0
;
i
<
NumPart
;
i
++
)
for
(
j
=
0
;
j
<
3
;
j
++
)
{
while
(
P
[
i
].
Pos
[
j
]
<
0
)
P
[
i
].
Pos
[
j
]
+=
boxsize
[
j
];
while
(
P
[
i
].
Pos
[
j
]
>=
boxsize
[
j
])
P
[
i
].
Pos
[
j
]
-=
boxsize
[
j
];
}
}
#endif
Event Timeline
Log In to Comment