testing/icgen/random_verschillende_resoluties_N-GenIC/gadget3_64/timestep.c
void find_dt_displacement_constraint(double hfac /*!<  should be  a^2*H(a)  */ )
{
  int i, type;
  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(P[i].Mass > 0)
	    {
	      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);

      sumup_large_ints(6, count, count_sum);

#ifdef SFR
      /* add star and gas particles together to treat them on equal footing, using the original gas particle
         spacing. */
      v_sum[0] += v_sum[4];
      count_sum[0] += count_sum[4];
      v_sum[4] = v_sum[0];
      count_sum[4] = count_sum[0];
#ifdef BLACK_HOLES
      v_sum[0] += v_sum[5];
      count_sum[0] += count_sum[5];
      v_sum[5] = v_sum[0];
      count_sum[5] = count_sum[0];
      min_mass[5] = min_mass[0];
#endif
#endif

      for(type = 0; type < 6; type++)
	{
	  if(count_sum[type] > 0)
	    {
	      if(type == 0 || (type == 4 && All.StarformationOn))
		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);

#ifdef BLACK_HOLES
	      if(type == 5)
		dmean =
		  pow(min_mass[type] / (All.OmegaBaryon * 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G)),
		      1.0 / 3);
#endif
	      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);


#ifdef NEUTRINOS
	      if(type != 2)	/* don't constrain the step to the neutrinos */
#endif
		if(dt < dt_displacement)
		  dt_displacement = dt;
	    }
	}

      if(ThisTask == 0)
	printf("displacement time constraint: %g  (%g)\n", dt_displacement, All.MaxSizeTimestep);
    }
}
/*! 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.
 */
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.
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  )
{
  int i, type;
  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(P[i].Mass > 0)
	    {
	      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);
      sumup_large_ints(6, count, count_sum);
#ifdef SFR
      v_sum[0] += v_sum[4];
      count_sum[0] += count_sum[4];
      v_sum[4] = v_sum[0];
      count_sum[4] = count_sum[0];
#ifdef BLACK_HOLES
      v_sum[0] += v_sum[5];
      count_sum[0] += count_sum[5];
      v_sum[5] = v_sum[0];
      count_sum[5] = count_sum[0];
      min_mass[5] = min_mass[0];
#endif
#endif
      for(type = 0; type < 6; type++)
	{
	  if(count_sum[type] > 0)
	    {
	      if(type == 0 || (type == 4 && All.StarformationOn))
		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);
#ifdef BLACK_HOLES
	      if(type == 5)
		dmean =
		  pow(min_mass[type] / (All.OmegaBaryon * 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G)),
		      1.0 / 3);
#endif
	      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);
#ifdef NEUTRINOS
	      if(type != 2)	
#endif
		if(dt < dt_displacement)
		  dt_displacement = dt;
	    }
	}
      if(ThisTask == 0)
	printf("displacement time constraint: %g  (%g)\n", dt_displacement, All.MaxSizeTimestep);
    }
}
void find_dt_displacement_constraint ( double hfac ) { int i , type ; 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 ( P [ i ] . Mass > 0 ) { 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 ) ; sumup_large_ints ( 6 , count , count_sum ) ; #ifdef SFR v_sum [ 0 ] += v_sum [ 4 ] ; count_sum [ 0 ] += count_sum [ 4 ] ; v_sum [ 4 ] = v_sum [ 0 ] ; count_sum [ 4 ] = count_sum [ 0 ] ; #ifdef BLACK_HOLES v_sum [ 0 ] += v_sum [ 5 ] ; count_sum [ 0 ] += count_sum [ 5 ] ; v_sum [ 5 ] = v_sum [ 0 ] ; count_sum [ 5 ] = count_sum [ 0 ] ; min_mass [ 5 ] = min_mass [ 0 ] ; #endif #endif for ( type = 0 ; type < 6 ; type ++ ) { if ( count_sum [ type ] > 0 ) { if ( type == 0 || ( type == 4 && All . StarformationOn ) ) 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 ) ; #ifdef BLACK_HOLES if ( type == 5 ) dmean = pow ( min_mass [ type ] / ( All . OmegaBaryon * 3 * All . Hubble * All . Hubble / ( 8 * M_PI * All . G ) ) , 1.0 / 3 ) ; #endif 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 ( " \n " , type , dmean , asmth , min_mass [ type ] , All . Time , sqrt ( v_sum [ type ] / count_sum [ type ] ) , dt ) ; #ifdef NEUTRINOS if ( type != 2 ) #endif if ( dt < dt_displacement ) dt_displacement = dt ; } } if ( ThisTask == 0 ) printf ( " \n " , dt_displacement , All . MaxSizeTimestep ) ; } }
This function computes an upper limit ('dt_displacement') to the global timestep of the system based on
  the rms velocities of particles.
This function computes an upper limit ( ' dt_displacement ' ) to the global timestep of the system based on the rms velocities of particles .
/*!<  should be  a^2*H(a)  */ /* add star and gas particles together to treat them on equal footing, using the original gas particle
         spacing. */ /* don't constrain the step to the neutrinos */