Homec4science
Diffusion Gear (master) Browse

Gear/src/timestep.c
b87f34a0ce7emaster

Last Change
Raw File

timestep.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 timestep.c
* \brief routines for 'kicking' particles in momentum space and assigning new timesteps
*/
static double fac1, fac2, fac3, hubble_a, atime, a3inv;
static double dt_displacement = 0;
/*! This function advances the system in momentum space, i.e. it does apply
* the 'kick' operation after the forces have been computed. Additionally, it
* assigns new timesteps to particles. At start-up, a half-timestep is
* carried out, as well as at the end of the simulation. In between, the
* half-step kick that ends the previous timestep and the half-step kick for
* the new timestep are combined into one operation.
*/
void advance_and_find_timesteps(void)
{
int i, j, no, ti_step, ti_min, tend, tstart;
double dt_entr, dt_entr2, dt_gravkick, dt_hydrokick, dt_gravkick2, dt_hydrokick2, t0, t1;
double minentropy, aphys;
FLOAT dv[3];
#ifdef COOLING
double t2,t3;
#endif
#ifdef FLEXSTEPS
int ti_grp;
#endif
#if defined(PSEUDOSYMMETRIC) && !defined(FLEXSTEPS)
double apred, prob;
int ti_step2;
#endif
#ifdef PMGRID
double dt_gravkickA, dt_gravkickB;
#endif
#ifdef MAKEGLASS
double disp, dispmax, globmax, dmean, fac, disp2sum, globdisp2sum;
#endif
t0 = second();
if(All.ComovingIntegrationOn)
{
fac1 = 1 / (All.Time * All.Time);
fac2 = 1 / pow(All.Time, 3 * GAMMA - 2);
fac3 = pow(All.Time, 3 * (1 - GAMMA) / 2.0);
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);
a3inv = 1 / (All.Time * All.Time * All.Time);
atime = All.Time;
}
else
fac1 = fac2 = fac3 = hubble_a = a3inv = atime = 1;
#ifdef COOLING_GRACKLE
double a_now;
if(All.ComovingIntegrationOn)
a_now = All.Time;
else
{
if (All.GrackleRedshift>=0)
/* compute scaling factor from the GrackleRedshift parameter */
a_now = 1. / (1. + All.GrackleRedshift);
else
/* compute scaling factor from cosmic time */
a_now = get_a_from_CosmicTime(All.Time);
}
if(ThisTask==0)
printf("Set Grackle scaling factor to = %g\n",a_now);
#endif
#ifdef NOPMSTEPADJUSTMENT
dt_displacement = All.MaxSizeTimestep;
#else
if(Flag_FullStep || dt_displacement == 0)
find_dt_displacement_constraint(hubble_a * atime * atime);
#endif
#ifdef PMGRID
if(All.ComovingIntegrationOn)
dt_gravkickB = get_gravkick_factor(All.PM_Ti_begstep, All.Ti_Current) -
get_gravkick_factor(All.PM_Ti_begstep, (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2);
else
dt_gravkickB = (All.Ti_Current - (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2) * All.Timebase_interval;
if(All.PM_Ti_endstep == All.Ti_Current) /* need to do long-range kick */
{
/* make sure that we reconstruct the domain/tree next time because we don't kick the tree nodes in this case */
All.NumForcesSinceLastDomainDecomp = 1 + All.TotNumPart * All.TreeDomainUpdateFrequency;
}
#endif
#ifdef MAKEGLASS
for(i = 0, dispmax = 0, disp2sum = 0; i < NumPart; i++)
{
for(j = 0; j < 3; j++)
{
P[i].GravPM[j] *= -1;
P[i].GravAccel[j] *= -1;
P[i].GravAccel[j] += P[i].GravPM[j];
P[i].GravPM[j] = 0;
}
disp = sqrt(P[i].GravAccel[0] * P[i].GravAccel[0] +
P[i].GravAccel[1] * P[i].GravAccel[1] + P[i].GravAccel[2] * P[i].GravAccel[2]);
disp *= 2.0 / (3 * All.Hubble * All.Hubble);
disp2sum += disp * disp;
if(disp > dispmax)
dispmax = disp;
}
MPI_Allreduce(&dispmax, &globmax, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
MPI_Allreduce(&disp2sum, &globdisp2sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
dmean = pow(P[0].Mass / (All.Omega0 * 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G)), 1.0 / 3);
if(globmax > dmean)
fac = dmean / globmax;
else
fac = 1.0;
if(ThisTask == 0)
{
printf("\nglass-making: dmean= %g global disp-maximum= %g rms= %g\n\n",
dmean, globmax, sqrt(globdisp2sum / All.TotNumPart));
fflush(stdout);
}
for(i = 0, dispmax = 0; i < NumPart; i++)
{
for(j = 0; j < 3; j++)
{
P[i].Vel[j] = 0;
P[i].Pos[j] += fac * P[i].GravAccel[j] * 2.0 / (3 * All.Hubble * All.Hubble);
P[i].GravAccel[j] = 0;
}
}
#endif
/* Now assign new timesteps and kick */
#ifdef FLEXSTEPS
if((All.Ti_Current % (4 * All.PresentMinStep)) == 0)
if(All.PresentMinStep < TIMEBASE)
All.PresentMinStep *= 2;
for(i = 0; i < NumPart; i++)
{
if(P[i].Ti_endstep == All.Ti_Current)
{
ti_step = get_timestep(i, &aphys, 0);
/* make it a power 2 subdivision */
ti_min = TIMEBASE;
while(ti_min > ti_step)
ti_min >>= 1;
ti_step = ti_min;
if(ti_step < All.PresentMinStep)
All.PresentMinStep = ti_step;
}
}
ti_step = All.PresentMinStep;
MPI_Allreduce(&ti_step, &All.PresentMinStep, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
if(dt_displacement < All.MaxSizeTimestep)
ti_step = (int) (dt_displacement / All.Timebase_interval);
else
ti_step = (int) (All.MaxSizeTimestep / All.Timebase_interval);
/* make it a power 2 subdivision */
ti_min = TIMEBASE;
while(ti_min > ti_step)
ti_min >>= 1;
All.PresentMaxStep = ti_min;
if(ThisTask == 0)
printf("Syn Range = %g PresentMinStep = %d PresentMaxStep = %d \n",
(double) All.PresentMaxStep / All.PresentMinStep, All.PresentMinStep, All.PresentMaxStep);
#endif
#ifdef SYNCHRONIZE_NGB_TIMESTEP
for(i = 0; i < NumPart; i++)
{
P[i].Old_Ti_begstep = P[i].Ti_begstep;
P[i].Old_Ti_endstep = P[i].Ti_endstep;
}
#endif
for(i = 0; i < NumPart; i++)
{
if(P[i].Ti_endstep == All.Ti_Current)
{
ti_step = get_timestep(i, &aphys, 0);
/* make it a power 2 subdivision */
ti_min = TIMEBASE;
while(ti_min > ti_step)
ti_min >>= 1;
ti_step = ti_min;
#ifdef FLEXSTEPS
ti_grp = P[i].FlexStepGrp % All.PresentMaxStep;
ti_grp = (ti_grp / All.PresentMinStep) * All.PresentMinStep;
ti_step = ((P[i].Ti_endstep + ti_grp + ti_step) / ti_step) * ti_step - (P[i].Ti_endstep + ti_grp);
#else
#ifdef PSEUDOSYMMETRIC
if(P[i].Type != 0)
{
if(P[i].Ti_endstep > P[i].Ti_begstep)
{
apred = aphys + ((aphys - P[i].AphysOld) / (P[i].Ti_endstep - P[i].Ti_begstep)) * ti_step;
if(fabs(apred - aphys) < 0.5 * aphys)
{
ti_step2 = get_timestep(i, &apred, -1);
ti_min = TIMEBASE;
while(ti_min > ti_step2)
ti_min >>= 1;
ti_step2 = ti_min;
if(ti_step2 < ti_step)
{
get_timestep(i, &apred, ti_step);
prob =
((apred - aphys) / (aphys - P[i].AphysOld) * (P[i].Ti_endstep -
P[i].Ti_begstep)) / ti_step;
if(prob < get_random_number(P[i].ID))
ti_step /= 2;
}
else if(ti_step2 > ti_step)
{
get_timestep(i, &apred, 2 * ti_step);
prob =
((apred - aphys) / (aphys - P[i].AphysOld) * (P[i].Ti_endstep -
P[i].Ti_begstep)) / ti_step;
if(prob < get_random_number(P[i].ID + 1))
ti_step *= 2;
}
}
}
P[i].AphysOld = aphys;
}
#endif
#ifdef SYNCHRONIZATION
if(ti_step > (P[i].Ti_endstep - P[i].Ti_begstep)) /* timestep wants to increase */
{
//if(((TIMEBASE - P[i].Ti_endstep) % ti_step) > 0)
// ti_step = P[i].Ti_endstep - P[i].Ti_begstep; /* leave at old step */
while(((TIMEBASE - P[i].Ti_endstep) % ti_step) > 0) /* yr : allow to increase */
ti_step = ti_step/2;
}
#endif
#endif /* end of FLEXSTEPS */
if(All.Ti_Current == TIMEBASE) /* we here finish the last timestep. */
ti_step = 0;
if((TIMEBASE - All.Ti_Current) < ti_step) /* check that we don't run beyond the end */
ti_step = TIMEBASE - All.Ti_Current;
#ifdef SYNCHRONIZE_NGB_TIMESTEP
/* Here, in order to performe the synchronization of the time steps
* for neighbor particles, we need to interupt the loop.
*/
P[i].Ti_step = ti_step; /* save estimated time step */
}
}
synchronize_ngb_timestep();
for(i = 0; i < NumPart; i++)
{
if(P[i].Old_Ti_endstep == All.Ti_Current) // here we use the old value, avoid problem due to the update of the timesteps
{
ti_step = P[i].Ti_step; /* recover from the estimated time step */
#endif
tstart = (P[i].Ti_begstep + P[i].Ti_endstep) / 2; /* midpoint of old step */
tend = P[i].Ti_endstep + ti_step / 2; /* midpoint of new step */
if(All.ComovingIntegrationOn)
{
dt_entr = (tend - tstart) * All.Timebase_interval;
dt_entr2 = (tend - P[i].Ti_endstep) * All.Timebase_interval;
dt_gravkick = get_gravkick_factor(tstart, tend);
dt_hydrokick = get_hydrokick_factor(tstart, tend);
dt_gravkick2 = get_gravkick_factor(P[i].Ti_endstep, tend);
dt_hydrokick2 = get_hydrokick_factor(P[i].Ti_endstep, tend);
}
else
{
dt_entr = dt_gravkick = dt_hydrokick = (tend - tstart) * All.Timebase_interval;
dt_gravkick2 = dt_hydrokick2 = dt_entr2 = (tend - P[i].Ti_endstep) * All.Timebase_interval;
}
P[i].Ti_begstep = P[i].Ti_endstep;
P[i].Ti_endstep = P[i].Ti_begstep + ti_step;
#ifdef CYLINDRICAL_SYMMETRY
double r,factor;
r = sqrt( P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2] );
factor = 1/(r*r) * (P[i].Pos[0]*P[i].GravAccel[0] + P[i].Pos[1]*P[i].GravAccel[1]);
P[i].GravAccel[0] = factor * P[i].Pos[0];
P[i].GravAccel[1] = factor * P[i].Pos[1];
#endif
#ifdef SIDM_ON
if(P[i].Type == SIDM_type)
P[i].Dt = ti_step * All.Timebase_interval;
#endif
/* do the kick */
for(j = 0; j < 3; j++)
{
dv[j] = 0.0;
#ifdef LIMIT_DVEL
if (fabs(P[i].GravAccel[j] * dt_gravkick)>LIMIT_DVEL)
{
#ifdef MULTIPHASE
printf("Warning(LIMIT_DVEL): ID=%d j=%d dv[j]=%g Phase=%d(setting GravAccel[j] to 0.0)\n",P[i].ID,j,P[i].GravAccel[j]*dt_hydrokick,SphP[i].Phase);
#else
printf("Warning(LIMIT_DVEL): ID=%d j=%d dv[j]=%g Phase=-(setting GravAccel[j] to 0.0)\n",P[i].ID,j,P[i].GravAccel[j]*dt_hydrokick);
#endif
P[i].GravAccel[j] = 0.0;
}
#endif
dv[j] += P[i].GravAccel[j] * dt_gravkick;
P[i].Vel[j] += P[i].GravAccel[j] * dt_gravkick;
}
if(P[i].Type == 0) /* SPH stuff */
{
for(j = 0; j < 3; j++)
{
#ifdef LIMIT_DVEL /* begin LIMIT_DVEL */
if (fabs(SphP[i].HydroAccel[j] * dt_hydrokick)>LIMIT_DVEL)
{
#ifdef MULTIPHASE
printf("Warning(LIMIT_DVEL): ID=%d j=%d dv[j]=%g Phase=%d(setting HydroAccel[j] to 0.0)\n",P[i].ID,j,SphP[i].HydroAccel[j] *dt_hydrokick,SphP[i].Phase);
#else
printf("Warning(LIMIT_DVEL): ID=%d j=%d dv[j]=%g Phase=-(setting HydroAccel[j] to 0.0)\n",P[i].ID,j,SphP[i].HydroAccel[j] *dt_hydrokick);
#endif
SphP[i].HydroAccel[j] = 0.0;
}
#endif /* end LIMIT_DVEL */
dv[j] += SphP[i].HydroAccel[j] * dt_hydrokick;
P[i].Vel[j] += SphP[i].HydroAccel[j] * dt_hydrokick;
SphP[i].VelPred[j] =
P[i].Vel[j] - dt_gravkick2 * P[i].GravAccel[j] - dt_hydrokick2 * SphP[i].HydroAccel[j];
#ifdef PMGRID
SphP[i].VelPred[j] += P[i].GravPM[j] * dt_gravkickB;
#endif
#ifdef AB_TURB
dv[j] += SphP[i].TurbAccel[j] * dt_hydrokick;
P[i].Vel[j] += SphP[i].TurbAccel[j] * dt_hydrokick;
SphP[i].VelPred[j] += - dt_hydrokick2 * SphP[i].TurbAccel[j];
#endif
#ifdef DISSIPATION_FORCES
dv[j] += SphP[i].DissipationForcesAccel[j] * dt_hydrokick;
P[i].Vel[j] += SphP[i].DissipationForcesAccel[j] * dt_hydrokick;
SphP[i].VelPred[j] += - dt_hydrokick2 * SphP[i].DissipationForcesAccel[j] ;
#endif
}
/***********************************************************/
/* compute spec energy lost/win by different other process */
/***********************************************************/
/***********************************************************/
/* compute entropy variation */
/***********************************************************/
/*******************************/
/* compute cooling */
/*******************************/
#ifdef COOLING
t2 = second();
#ifdef COOLING_GRACKLE
CoolingForOne_GRACKLE(i,tstart,tend,dt_entr,dt_entr2,ti_step / 2 * All.Timebase_interval,a_now,a3inv);
#else
CoolingForOne(i,tstart,tend,ti_step,dt_entr2,a3inv,hubble_a);
#endif
t3 = second();
All.CPU_Cooling += timediff(t2, t3);
#else
/* In case of cooling, we prevent that the entropy (and
hence temperature decreases by more than a factor 0.5 */
if(SphP[i].DtEntropy * dt_entr > -0.5 * SphP[i].Entropy)
SphP[i].Entropy += SphP[i].DtEntropy * dt_entr;
else
SphP[i].Entropy *= 0.5;
#ifdef MULTIPHASE
if (SphP[i].Phase==GAS_SPH)
{
#endif
if(All.MinEgySpec)
{
#ifdef DENSITY_INDEPENDENT_SPH
minentropy = All.MinEgySpec * GAMMA_MINUS1 / pow(SphP[i].EgyWtDensity * a3inv, GAMMA_MINUS1);
#else
minentropy = All.MinEgySpec * GAMMA_MINUS1 / pow(SphP[i].Density * a3inv, GAMMA_MINUS1);
#endif
if(SphP[i].Entropy < minentropy)
{
SphP[i].Entropy = minentropy;
SphP[i].DtEntropy = 0;
}
}
#ifdef MULTIPHASE
}
#endif
#endif /* COOLING */
#ifndef COOLING
/* In case the timestep increases in the new step, we
make sure that we do not 'overcool' when deriving
predicted temperatures. The maximum timespan over
which prediction can occur is ti_step/2, i.e. from
the middle to the end of the current step */
//dt_entr = ti_step / 2 * All.Timebase_interval;
dt_entr = imax(ti_step / 2,1) * All.Timebase_interval; /* yr : prevent dt_entr to be zero if ti_step=1 */
if(SphP[i].Entropy + SphP[i].DtEntropy * dt_entr < 0.5 * SphP[i].Entropy)
SphP[i].DtEntropy = -0.5 * SphP[i].Entropy / dt_entr;
#ifdef ENTROPYPRED
/* now, we correct the predicted Entropy */
SphP[i].EntropyPred = SphP[i].Entropy - dt_entr2 * SphP[i].DtEntropy ;
#ifdef COOLING
// if(All.MinEgySpec)
// minentropy = All.MinEgySpec * GAMMA_MINUS1 / pow(SphP[i].Density * a3inv, GAMMA_MINUS1);
// else
// minentropy=0;
//
// //dt_entr = imax(ti_step,1) * All.Timebase_interval;
// //if (SphP[i].EntropyPred + SphP[i].DtEntropy * dt_entr < minentropy) /* if during the next step prediction entropy may be less than zero */
// // SphP[i].DtEntropy = -(SphP[i].EntropyPred-minentropy) / dt_entr; /* modify SphP[i].DtEntropy in order to avoid problems */
//
// dt_entr = imax(ti_step / 2,1) * All.Timebase_interval;
// if(SphP[i].Entropy + SphP[i].DtEntropy * dt_entr < 0.5 * SphP[i].Entropy)
// {
// printf("(%d) particle id=%d reduces its DtEntropy (Entropy=%g DtEntropy * dt_entr=%g)\n",NTask,P[i].ID,SphP[i].Entropy,SphP[i].DtEntropy * dt_entr);
// SphP[i].DtEntropy = -0.5 * SphP[i].Entropy / dt_entr;
// }
#endif /* COOLING */
#ifdef CHECK_ENTROPY_SIGN
if ((SphP[i].EntropyPred < 0))
{
printf("\ntask=%d: EntropyPred less than zero in advance_and_find_timesteps !\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(1010101000);
}
#endif
#endif /* ENTROPYPRED */
#endif /* no COOLING */
#ifdef DISSIPATION_FORCES
dt_entr = imax(ti_step / 2,1) * All.Timebase_interval; /* yr : prevent dt_entr to be zero if ti_step=1 */
SphP[i].EnergyDissipationForces += SphP[i].DtEnergyDissipationForces * dt_entr;
#endif
}
/* if tree is not going to be reconstructed, kick parent nodes dynamically.
*/
if(All.NumForcesSinceLastDomainDecomp < All.TotNumPart * All.TreeDomainUpdateFrequency)
{
no = Father[i];
while(no >= 0)
{
for(j = 0; j < 3; j++)
Extnodes[no].vs[j] += dv[j] * P[i].Mass / Nodes[no].u.d.mass;
no = Nodes[no].u.d.father;
}
}
}
}
#ifdef PMGRID
if(All.PM_Ti_endstep == All.Ti_Current) /* need to do long-range kick */
{
ti_step = TIMEBASE;
while(ti_step > (dt_displacement / All.Timebase_interval))
ti_step >>= 1;
if(ti_step > (All.PM_Ti_endstep - All.PM_Ti_begstep)) /* PM-timestep wants to increase */
{
/* we only increase if an integer number of steps will bring us to the end */
if(((TIMEBASE - All.PM_Ti_endstep) % ti_step) > 0)
ti_step = All.PM_Ti_endstep - All.PM_Ti_begstep; /* leave at old step */
}
if(All.Ti_Current == TIMEBASE) /* we here finish the last timestep. */
ti_step = 0;
tstart = (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2;
tend = All.PM_Ti_endstep + ti_step / 2;
if(All.ComovingIntegrationOn)
dt_gravkick = get_gravkick_factor(tstart, tend);
else
dt_gravkick = (tend - tstart) * All.Timebase_interval;
All.PM_Ti_begstep = All.PM_Ti_endstep;
All.PM_Ti_endstep = All.PM_Ti_begstep + ti_step;
if(All.ComovingIntegrationOn)
dt_gravkickB = -get_gravkick_factor(All.PM_Ti_begstep, (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2);
else
dt_gravkickB =
-((All.PM_Ti_begstep + All.PM_Ti_endstep) / 2 - All.PM_Ti_begstep) * All.Timebase_interval;
for(i = 0; i < NumPart; i++)
{
for(j = 0; j < 3; j++) /* do the kick */
P[i].Vel[j] += P[i].GravPM[j] * dt_gravkick;
if(P[i].Type == 0)
{
if(All.ComovingIntegrationOn)
{
dt_gravkickA = get_gravkick_factor(P[i].Ti_begstep, All.Ti_Current) -
get_gravkick_factor(P[i].Ti_begstep, (P[i].Ti_begstep + P[i].Ti_endstep) / 2);
dt_hydrokick = get_hydrokick_factor(P[i].Ti_begstep, All.Ti_Current) -
get_hydrokick_factor(P[i].Ti_begstep, (P[i].Ti_begstep + P[i].Ti_endstep) / 2);
}
else
dt_gravkickA = dt_hydrokick =
(All.Ti_Current - (P[i].Ti_begstep + P[i].Ti_endstep) / 2) * All.Timebase_interval;
for(j = 0; j < 3; j++)
SphP[i].VelPred[j] = P[i].Vel[j]
+ P[i].GravAccel[j] * dt_gravkickA
+ SphP[i].HydroAccel[j] * dt_hydrokick
+ P[i].GravPM[j] * dt_gravkickB;
}
}
}
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
chimie_apply_thermal_feedback();
#endif
t1 = second();
All.CPU_TimeLine += timediff(t0, t1);
#ifdef DETAILED_CPU
All.CPU_Leapfrog += timediff(t0, t1);
#endif
#ifdef COOLING
//All.CPU_TimeLine -= All.CPU_Cooling;
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
if(SetMinTimeStepForActives)
SetMinTimeStepForActives=0;
#endif
#ifdef SYNCHRONIZE_NGB_TIMESTEP
#ifdef OUTPUTOPTVAR1
for(i = 0; i < N_gas; i++)
SphP[i].OptVar1 = (float) (P[i].Ti_endstep - P[i].Ti_begstep);
#endif
#endif
}
/*! This function normally (for flag==0) returns the maximum allowed timestep
* of a particle, expressed in terms of the integer mapping that is used to
* represent the total simulated timespan. The physical acceleration is
* returned in `aphys'. The latter is used in conjunction with the
* PSEUDOSYMMETRIC integration option, which also makes of the second
* function of get_timestep. When it is called with a finite timestep for
* flag, it returns the physical acceleration that would lead to this
* timestep, assuming timestep criterion 0.
*/
int get_timestep(int p, /*!< particle index */
double *aphys, /*!< acceleration (physical units) */
int flag /*!< either 0 for normal operation, or finite timestep to get corresponding
aphys */ )
{
double ax, ay, az, ac, csnd;
double dt = 0, dt_courant = 0, dt_accel;
double maxsignalvel;
int ti_step;
#ifdef CONDUCTION
double dt_cond;
#endif
if(flag == 0)
{
ax = fac1 * P[p].GravAccel[0];
ay = fac1 * P[p].GravAccel[1];
az = fac1 * P[p].GravAccel[2];
#ifdef PMGRID
ax += fac1 * P[p].GravPM[0];
ay += fac1 * P[p].GravPM[1];
az += fac1 * P[p].GravPM[2];
#endif
if(P[p].Type == 0)
{
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
if((SphP[p].DeltaEgySpec==-1) || (SphP[p].DeltaEgySpec>0))
{
ax += fac2 * SphP[p].FeedbackUpdatedAccel[0];
ay += fac2 * SphP[p].FeedbackUpdatedAccel[1];
az += fac2 * SphP[p].FeedbackUpdatedAccel[2];
if(SphP[p].DeltaEgySpec==-1)
SphP[p].DeltaEgySpec=0; /* unflag */
}
else
{
ax += fac2 * SphP[p].HydroAccel[0];
ay += fac2 * SphP[p].HydroAccel[1];
az += fac2 * SphP[p].HydroAccel[2];
}
#else
ax += fac2 * SphP[p].HydroAccel[0];
ay += fac2 * SphP[p].HydroAccel[1];
az += fac2 * SphP[p].HydroAccel[2];
#endif
#ifdef DISSIPATION_FORCES
ax += fac2 * SphP[p].DissipationForcesAccel[0];
ay += fac2 * SphP[p].DissipationForcesAccel[1];
az += fac2 * SphP[p].DissipationForcesAccel[2];
#endif
#ifdef AB_TURB
ax += fac2 * SphP[p].TurbAccel[0];
ay += fac2 * SphP[p].TurbAccel[1];
az += fac2 * SphP[p].TurbAccel[2];
#endif
}
ac = sqrt(ax * ax + ay * ay + az * az); /* this is now the physical acceleration */
*aphys = ac;
}
else
ac = *aphys;
if(ac == 0)
ac = 1.0e-30;
switch (All.TypeOfTimestepCriterion)
{
case 0:
if(flag > 0)
{
dt = flag * All.Timebase_interval;
dt /= hubble_a; /* convert dloga to physical timestep */
ac = 2 * All.ErrTolIntAccuracy * atime * All.SofteningTable[P[p].Type] / (dt * dt);
*aphys = ac;
return flag;
}
#ifdef IMPROVED_TIMESTEP_CRITERION_FORGAS
if(P[p].Type == 0)
dt = dt_accel = sqrt(2 * All.ErrTolIntAccuracy * atime * dmin(SphP[p].Hsml, All.SofteningTable[P[p].Type]) / ac);
else
dt = dt_accel = sqrt(2 * All.ErrTolIntAccuracy * atime * All.SofteningTable[P[p].Type] / ac);
#else
dt = dt_accel = sqrt(2 * All.ErrTolIntAccuracy * atime * All.SofteningTable[P[p].Type] / ac);
#endif
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(P[p].Type == 0)
dt = dt_accel = sqrt(2 * All.ErrTolIntAccuracy * atime * SphP[p].Hsml / 2.8 / ac);
#endif
break;
default:
endrun(888);
break;
}
if(P[p].Type == 0)
{
#ifdef DENSITY_INDEPENDENT_SPH
csnd = sqrt(GAMMA * SphP[p].Pressure / SphP[p].EgyWtDensity);
#else
csnd = sqrt(GAMMA * SphP[p].Pressure / SphP[p].Density);
#endif
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
maxsignalvel = dmax(SphP[p].MaxSignalVel,SphP[p].MaxSignalVelFeedbackUpdated);
#else
maxsignalvel = SphP[p].MaxSignalVel;
#endif
if(All.ComovingIntegrationOn)
dt_courant = 2 * All.CourantFac * All.Time * SphP[p].Hsml / (fac3 * maxsignalvel);
else
dt_courant = 2 * All.CourantFac * SphP[p].Hsml / maxsignalvel;
if(dt_courant < dt)
#ifndef MULTIPHASE
dt = dt_courant;
#else
{
if (maxsignalvel != 0);
dt = dt_courant;
}
#endif
// #ifdef CHIMIE_THERMAL_FEEDBACK
//
// float f;
// double EgySpec,NewEgySpec;
//
// if (SphP[p].DeltaEgySpec > 0)
// {
//
// /* spec energy at current step */
// EgySpec = SphP[p].EntropyPred / GAMMA_MINUS1 * pow(SphP[p].Density*a3inv, GAMMA_MINUS1);
//
// /* new egyspec */
// NewEgySpec = EgySpec + SphP[p].DeltaEgySpec;
//
// f = NewEgySpec/EgySpec;
//
// //if (f>1)
// // dt = dt / f;
// }
//
// #endif
#ifdef CHIMIE_KINETIC_FEEDBACK
double dt_kinetic_feedback;
double SupernovaKieticFeedbackIntAccuracy=0.1;
dt_kinetic_feedback = SupernovaKieticFeedbackIntAccuracy * All.SofteningTable[P[p].Type] / All.ChimieWindSpeed;
if(dt_kinetic_feedback < dt)
dt = dt_kinetic_feedback;
#endif
#ifdef FEEDBACK_WIND
double dt_feedback_wind;
double vwind;
vwind = sqrt( SphP[p].FeedbackWindVel[0]*SphP[p].FeedbackWindVel[0] + SphP[p].FeedbackWindVel[1]*SphP[p].FeedbackWindVel[1] + SphP[p].FeedbackWindVel[2]*SphP[p].FeedbackWindVel[2] );
if (vwind > 0)
{
dt_feedback_wind = All.SupernovaWindIntAccuracy * All.SofteningTable[P[p].Type] / vwind;
SphP[p].FeedbackWindVel[0]=0;
SphP[p].FeedbackWindVel[1]=0;
SphP[p].FeedbackWindVel[2]=0;
if(dt_feedback_wind < dt)
dt = dt_feedback_wind;
}
#endif
}
#ifdef CHIMIE
int m;
double dt_chimie;
if(P[p].Type == ST)
{
m = P[p].StPIdx;
if (StP[m].Flag)
{
if(All.ComovingIntegrationOn)
dt_chimie = All.ChimieMaxSizeTimestep * All.Time * hubble_a;
else
dt_chimie = All.ChimieMaxSizeTimestep;
if(dt_chimie < dt)
dt = dt_chimie;
}
}
#endif
/* convert the physical timestep to dloga if needed. Note: If comoving integration has not been selected,
hubble_a=1.
*/
dt *= hubble_a;
if(dt >= All.MaxSizeTimestep)
dt = All.MaxSizeTimestep;
if(dt >= dt_displacement)
dt = dt_displacement;
if(dt < All.MinSizeTimestep)
{
#ifndef NOSTOP_WHEN_BELOW_MINTIMESTEP
printf("warning: Timestep wants to be below the limit `MinSizeTimestep'\n");
if(P[p].Type == 0)
{
printf
("Part-ID=%d dt=%g dtc=%g ac=%g xyz=(%g|%g|%g) hsml=%g maxsignalvel=%g dt0=%g eps=%g\n",
(int) P[p].ID, dt, dt_courant * hubble_a, ac, P[p].Pos[0], P[p].Pos[1], P[p].Pos[2],
SphP[p].Hsml, maxsignalvel,
sqrt(2 * All.ErrTolIntAccuracy * atime * All.SofteningTable[P[p].Type] / ac) * hubble_a,
All.SofteningTable[P[p].Type]);
}
else
{
printf("Part-ID=%d dt=%g ac=%g xyz=(%g|%g|%g)\n", (int) P[p].ID, dt, ac, P[p].Pos[0], P[p].Pos[1],
P[p].Pos[2]);
}
fflush(stdout);
endrun(888);
#endif
dt = All.MinSizeTimestep;
}
ti_step = dt / All.Timebase_interval;
#ifdef CHIMIE_KINETIC_FEEDBACK
//if(SetMinTimeStepForActives)
// ti_step=1;
#endif
if(!(ti_step > 0 && ti_step < TIMEBASE))
{
printf("\nError: A timestep of size zero was assigned on the integer timeline!\n"
"We better stop.\n"
"Task=%d Part-ID=%d dt=%g tibase=%g ti_step=%d ac=%g xyz=(%g|%g|%g) tree=(%g|%g%g)\n\n",
ThisTask, (int) P[p].ID, dt, All.Timebase_interval, ti_step, ac,
P[p].Pos[0], P[p].Pos[1], P[p].Pos[2], P[p].GravAccel[0], P[p].GravAccel[1], P[p].GravAccel[2]);
#ifdef PMGRID
printf("pm_force=(%g|%g|%g)\n", P[p].GravPM[0], P[p].GravPM[1], P[p].GravPM[2]);
#endif
if(P[p].Type == 0)
printf("hydro-frc=(%g|%g|%g)\n", SphP[p].HydroAccel[0], SphP[p].HydroAccel[1], SphP[p].HydroAccel[2]);
#ifdef FEEDBACK_WIND
if(P[p].Type == 0)
printf("feedback-vel=(%g|%g|%g)\n", SphP[p].FeedbackWindVel[0], SphP[p].FeedbackWindVel[1], SphP[p].FeedbackWindVel[2]);
#endif
fflush(stdout);
endrun(818);
}
return ti_step;
}
/*! This function computes an upper limit ('dt_displacement') to the global
* timestep of the system based on the rms velocities of particles. For
* cosmological simulations, the criterion used is that the rms displacement
* should be at most a fraction MaxRMSDisplacementFac of the mean particle
* separation. Note that the latter is estimated using the assigned particle
* masses, separately for each particle type. If comoving integration is not
* used, the function imposes no constraint on the timestep.
*/
void find_dt_displacement_constraint(double hfac /*!< should be a^2*H(a) */ )
{
int i, j, type, *temp;
int count[6];
long long count_sum[6];
double v[6], v_sum[6], mim[6], min_mass[6];
double dt, dmean, asmth = 0;
dt_displacement = All.MaxSizeTimestep;
if(All.ComovingIntegrationOn)
{
for(type = 0; type < 6; type++)
{
count[type] = 0;
v[type] = 0;
mim[type] = 1.0e30;
}
for(i = 0; i < NumPart; i++)
{
v[P[i].Type] += 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(mim[P[i].Type] > P[i].Mass)
mim[P[i].Type] = P[i].Mass;
count[P[i].Type]++;
}
MPI_Allreduce(v, v_sum, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(mim, min_mass, 6, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
temp = malloc(NTask * 6 * sizeof(int));
MPI_Allgather(count, 6, MPI_INT, temp, 6, MPI_INT, MPI_COMM_WORLD);
for(i = 0; i < 6; i++)
{
count_sum[i] = 0;
for(j = 0; j < NTask; j++)
count_sum[i] += temp[j * 6 + i];
}
free(temp);
for(type = 0; type < 6; type++)
{
if(count_sum[type] > 0)
{
if(type == 0)
dmean =
pow(min_mass[type] / (All.OmegaBaryon * 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G)),
1.0 / 3);
else
dmean =
pow(min_mass[type] /
((All.Omega0 - All.OmegaBaryon) * 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G)),
1.0 / 3);
dt = All.MaxRMSDisplacementFac * hfac * dmean / sqrt(v_sum[type] / count_sum[type]);
#ifdef PMGRID
asmth = All.Asmth[0];
#ifdef PLACEHIGHRESREGION
if(((1 << type) & (PLACEHIGHRESREGION)))
asmth = All.Asmth[1];
#endif
if(asmth < dmean)
dt = All.MaxRMSDisplacementFac * hfac * asmth / sqrt(v_sum[type] / count_sum[type]);
#endif
if(ThisTask == 0)
printf("type=%d dmean=%g asmth=%g minmass=%g a=%g sqrt(<p^2>)=%g dlogmax=%g\n",
type, dmean, asmth, min_mass[type], All.Time, sqrt(v_sum[type] / count_sum[type]), dt);
if(dt < dt_displacement)
dt_displacement = dt;
}
}
if(ThisTask == 0)
printf("displacement time constraint: %g (%g)\n", dt_displacement, All.MaxSizeTimestep);
}
}
#ifdef SYNCHRONIZE_NGB_TIMESTEP
#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
static int NDone;
static long long NTotDone;
/*! This function share the timesteps between particles
* according the the Saitoh and Makino rule.
*/
void synchronize_ngb_timestep()
{
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 t0, t1;
double timecomp = 0, timeimbalance = 0, timecommsumm = 0;
MPI_Status status;
int CptLimit = 0;
int shrinkcount = 0, shrinktot = 0;
int tstart,tend;
double dt_entr;
double dt_gravkick;
double dt_hydrokick;
int counter;
#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
/* `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
#ifdef MULTIPHASE
if(SphP[n].Phase == GAS_SPH)
#endif
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);
if (ThisTask==0)
printf("start synchronize ngb timestep\n");
NTotDone=1;
//NTotDone=0; /* if we want to disable the multi loop */
/* loop for time-step limiter */
while(NTotDone > 0)
{
NDone = 0;
i = 0; /* first particle for this task */
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 */
t0 = second();
for(nexport = 0, ndone = 0; i < N_gas && nexport < All.BunchSizeSynchronizeNgBTimestep - 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
{
{
ndone++;
for(j = 0; j < NTask; j++)
Exportflag[j] = 0;
NDone += synchronize_ngb_timestep_evaluate(i, 0);
for(j = 0; j < NTask; j++)
{
if(Exportflag[j])
{
SynchroinzeNgbTimestepDataIn[nexport].Pos[0] = P[i].Pos[0];
SynchroinzeNgbTimestepDataIn[nexport].Pos[1] = P[i].Pos[1];
SynchroinzeNgbTimestepDataIn[nexport].Pos[2] = P[i].Pos[2];
SynchroinzeNgbTimestepDataIn[nexport].Hsml = SphP[i].Hsml;
SynchroinzeNgbTimestepDataIn[nexport].Ti_step = P[i].Ti_step;
SynchroinzeNgbTimestepDataIn[nexport].Ti_endstep = P[i].Ti_endstep;
SynchroinzeNgbTimestepDataIn[nexport].Index = i;
SynchroinzeNgbTimestepDataIn[nexport].Task = j;
#ifdef MULTIPHASE
SynchroinzeNgbTimestepDataIn[nexport].Phase = SphP[i].Phase;
#endif
nexport++;
nsend_local[j]++;
}
}
}
}
t1 = second();
timecomp += timediff(t0, t1);
qsort(SynchroinzeNgbTimestepDataIn, nexport, sizeof(struct SynchroinzeNgbTimestepdata_in), synchronize_ngb_timestep_compare_key);
for(j = 1, noffset[0] = 0; j < NTask; j++)
noffset[j] = noffset[j - 1] + nsend_local[j - 1];
t0 = second();
MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);
t1 = second();
timeimbalance += timediff(t0, t1);
/* now do the particles that need to be exported */
for(level = 1; level < (1 << PTask); level++)
{
t0 = 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.BunchSizeSynchronizeNgBTimestep)
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(&SynchroinzeNgbTimestepDataIn[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct SynchroinzeNgbTimestepdata_in), MPI_BYTE,
recvTask, 0,
&SynchroinzeNgbTimestepDataGet[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct SynchroinzeNgbTimestepdata_in),
MPI_BYTE, recvTask, 0, MPI_COMM_WORLD, &status);
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
t1 = second();
timecommsumm += timediff(t0, t1);
t0 = second();
for(j = 0; j < nbuffer[ThisTask]; j++)
NDone += synchronize_ngb_timestep_evaluate(j, 1);
t1 = second();
timecomp += timediff(t0, t1);
/* do a block to explicitly measure imbalance */
t0 = second();
MPI_Barrier(MPI_COMM_WORLD);
t1 = second();
timeimbalance += timediff(t0, t1);
level = ngrp - 1;
}
t0 = second();
MPI_Allgather(&ndone, 1, MPI_INT, ndonelist, 1, MPI_INT, MPI_COMM_WORLD);
for(j = 0; j < NTask; j++)
ntotleft -= ndonelist[j];
t1 = second();
timeimbalance += timediff(t0, t1);
}
t0 = second();
numlist = (int*)malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&NDone, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, NTotDone = 0; i < NTask; i++)
NTotDone += numlist[i];
free(numlist);
t1 = second();
timeimbalance += timediff(t0, t1);
if(ThisTask == 0)
{
fprintf(stdout," %3d) ---> number of timestep shrinked gas neighbors: %6lld \n", CptLimit++, NTotDone);
fflush(stdout);
}
}
free(ndonelist);
free(nsend);
free(nsend_local);
free(nbuffer);
free(noffset);
/* do final operations on results */
counter=0;
for(i = 0; i < N_gas; i++)
#ifdef SFR
if((P[i].Type == 0))
#endif
{
if (P[i].Old_Ti_endstep != All.Ti_Current) /* the particle is inactive */
{
if ( (P[i].Old_Ti_endstep != P[i].Ti_endstep) || (P[i].Old_Ti_begstep != P[i].Ti_begstep) ) /* its timestep has been updated */
{
//printf("---------> %d %d %d %d\n",P[i].Old_Ti_endstep,P[i].Ti_endstep,P[i].Old_Ti_begstep,P[i].Ti_begstep);
/* need to extrapolate mid-step quantities */
counter++;
/* old mid step */
tstart = (P[i].Old_Ti_begstep + P[i].Old_Ti_endstep) / 2; /* midpoint of old step */
tend = (P[i].Ti_begstep + P[i].Ti_endstep) / 2; /* midpoint of new step */
/* now, do the kick */
kickback(i,tstart,tend);
}
}
}
if (counter!=0)
printf("(%d) %d passive particles have been updated \n",ThisTask,counter);
if (ThisTask==0)
printf("synchronize ngb timestep done.\n");
}
int synchronize_ngb_timestep_evaluate(int target, int mode)
{
int j, k, n, startnode, numngb_inbox;
double h, h2;
double r2, dx, dy, dz;
int phase=0;
FLOAT *pos;
int ti_step_i;
int endstep_i;
int CptShrink = 0;
if(mode == 0)
{
pos = P[target].Pos;
h = SphP[target].Hsml;
ti_step_i = P[target].Ti_step;
endstep_i = P[target].Ti_endstep; /* !!! */
#ifdef MULTIPHASE
phase = SphP[target].Phase;
#endif
}
else
{
pos = SynchroinzeNgbTimestepDataGet[target].Pos;
h = SynchroinzeNgbTimestepDataGet[target].Hsml;
ti_step_i = SynchroinzeNgbTimestepDataGet[target].Ti_step;
endstep_i = SynchroinzeNgbTimestepDataGet[target].Ti_endstep; /* !!! */
#ifdef MULTIPHASE
phase = SynchroinzeNgbTimestepDataGet[target].Phase;
#endif
}
h2 = h * h;
startnode = All.MaxPart;
do
{
numngb_inbox = ngb_treefind_variable(&pos[0], h, phase, &startnode);
for(n = 0; n < numngb_inbox; n++)
{
j = Ngblist[n];
dx = P[j].Pos[0] - pos[0];
dy = P[j].Pos[1] - pos[1];
dz = P[j].Pos[2] - 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)
{
if( P[j].Ti_endstep == All.Ti_Current ) /* the particle is active */
{
if(P[j].Ti_step > All.NgbFactorTimestep*ti_step_i )
{
CptShrink++;
P[j].Ti_step = All.NgbFactorTimestep*ti_step_i;
}
}
else /* the particle is not active */
{
if( P[j].Ti_endstep > All.Ti_Current + All.NgbFactorTimestep*ti_step_i )
{
CptShrink++;
P[j].Old_Ti_begstep = P[j].Ti_begstep;
P[j].Old_Ti_endstep = P[j].Ti_endstep;
P[j].Ti_step = All.NgbFactorTimestep*ti_step_i;
#ifdef SYNCHRONIZATION
/* find Ti_endstep*/
for(k = 0; P[j].Ti_begstep + k * P[j].Ti_step <= All.Ti_Current; k++);
P[j].Ti_endstep = P[j].Ti_begstep + k * P[j].Ti_step ;
P[j].Ti_begstep = P[j].Ti_endstep - P[j].Ti_step;
#else
P[j].Ti_endstep = All.Ti_Current + P[j].Ti_step;
P[j].Ti_begstep = P[j].Ti_endstep - P[j].Ti_step;
#endif
}
}
}
}
}
while(startnode >= 0);
return CptShrink;
}
/*! 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 synchronize_ngb_timestep_compare_key(const void *a, const void *b)
{
if(((struct SynchroinzeNgbTimestepdata_in *) a)->Task < (((struct SynchroinzeNgbTimestepdata_in *) b)->Task))
return -1;
if(((struct SynchroinzeNgbTimestepdata_in *) a)->Task > (((struct SynchroinzeNgbTimestepdata_in *) b)->Task))
return +1;
return 0;
}
#endif
c4science · Help