timestep.c
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, Sep 27, 17:01

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
	a_now = 1. / (1. + All.GrackleRedshift);

	printf("Grackle scaling factor = %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/(rr) (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

	/* 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,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;
	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)
	{

	#ifndef TIMESTEP_UPDATE_FOR_FEEDBACK
	ax += fac2 * SphP[p].HydroAccel[0];
	ay += fac2 * SphP[p].HydroAccel[1];
	az += fac2 * SphP[p].HydroAccel[2];
	#else
	#ifdef 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

	#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(All.ComovingIntegrationOn)
	dt_courant = 2 * All.CourantFac * All.Time * SphP[p].Hsml / (fac3 * SphP[p].MaxSignalVel);
	else
	dt_courant = 2 * All.CourantFac * SphP[p].Hsml / SphP[p].MaxSignalVel;

	if(dt_courant < dt)
	#ifndef MULTIPHASE
	dt = dt_courant;
	#else
	{
	if (SphP[p].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)
	{
	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, SphP[p].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^2H(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

timestep.cNo OneTemporaryActions

File Metadata

timestep.cView Options

Event Timeline

timestep.c
No OneTemporary
Actions

timestep.c
View Options