Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77547941
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
Thu, Aug 15, 02:58
Size
4 KB
Mime Type
text/x-c
Expires
Sat, Aug 17, 02:58 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
19886183
Attached To
R195 SCITAS application benchmarks
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;
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
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;
SphP[i].Pressure = (SphP[i].Entropy + SphP[i].DtEntropy * dt_entr) * pow(SphP[i].Density, GAMMA);
#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("\ntask=%d: entropy 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(333001);
}
#endif
#ifdef NO_NEGATIVE_PRESSURE
if (SphP[i].Pressure<0)
{
printf("\ntask=%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(333002);
}
#endif
/***********************************************************/
/* compute art visc coeff */
/***********************************************************/
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_ROSS) || defined(ART_VISCO_CD)
move_art_visc(i,dt_drift);
#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