testing/icgen/random_verschillende_resoluties_N-GenIC/gadget3_64/kinfb_mhm.c
void kinetic_evaluate(int target, int mode)
{
  int j, n;
  int startnode, numngb, numngb_inbox;
  double h, h2, hinv, hinv3;
  double rho, weight, wk, energy, delta_v;
  double dx, dy, dz, r, r2, u;
  double dvx, dvy, dvz;
  MyFloat *pos;

  if(mode == 0)
    {
      pos = P[target].Pos;
      h = PPP[target].Hsml;
      rho = SphP[target].Density;
      energy = SphP[target].FeedbackEnergy;
    }
  else
    {
      pos = KinDataGet[target].Pos;
      h = KinDataGet[target].Hsml;
      rho = KinDataGet[target].Density;
      energy = KinDataGet[target].Energy;
    }


  h2 = h * h;
  hinv = 1.0 / h;
#ifndef  TWODIMS
  hinv3 = hinv * hinv * hinv;
#else
  hinv3 = hinv * hinv / boxSize_Z;
#endif


  startnode = All.MaxPart;
  numngb = 0;
  do
    {
      numngb_inbox = ngb_treefind_variable(&pos[0], h, &startnode);

      for(n = 0; n < numngb_inbox; n++)
	{
	  j = Ngblist[n];

	  dx = pos[0] - P[j].Pos[0];
	  dy = pos[1] - P[j].Pos[1];
	  dz = pos[2] - P[j].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)
	    {
	      r = sqrt(r2);

	      if(r > 0)
		{

		  u = r * hinv;

		  if(u < 0.5)
		    wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
		  else
		    wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);

		  weight = wk * P[j].Mass / rho;

		  delta_v = sqrt(2 * energy * weight / P[j].Mass);

		  dvx = delta_v * (-dx) / r;
		  dvy = delta_v * (-dy) / r;
		  dvz = delta_v * (-dz) / r;

		  P[j].Vel[0] += dvx;
		  P[j].Vel[1] += dvy;
		  P[j].Vel[2] += dvz;

		  SphP[j].VelPred[0] += dvx;
		  SphP[j].VelPred[1] += dvy;
		  SphP[j].VelPred[2] += dvz;
		}
	    }
	}
    }
  while(startnode >= 0);
}
/*! This function represents the core of the SPH density computation. The
 *  target particle may either be local, or reside in the communication
 *  buffer.
 */
This function represents the core of the SPH density computation. The
target particle may either be local, or reside in the communication
buffer.
This function represents the core of the SPH density computation . The target particle may either be local or reside in the communication buffer .
void kinetic_evaluate(int target, int mode)
{
  int j, n;
  int startnode, numngb, numngb_inbox;
  double h, h2, hinv, hinv3;
  double rho, weight, wk, energy, delta_v;
  double dx, dy, dz, r, r2, u;
  double dvx, dvy, dvz;
  MyFloat *pos;
  if(mode == 0)
    {
      pos = P[target].Pos;
      h = PPP[target].Hsml;
      rho = SphP[target].Density;
      energy = SphP[target].FeedbackEnergy;
    }
  else
    {
      pos = KinDataGet[target].Pos;
      h = KinDataGet[target].Hsml;
      rho = KinDataGet[target].Density;
      energy = KinDataGet[target].Energy;
    }
  h2 = h * h;
  hinv = 1.0 / h;
#ifndef  TWODIMS
  hinv3 = hinv * hinv * hinv;
#else
  hinv3 = hinv * hinv / boxSize_Z;
#endif
  startnode = All.MaxPart;
  numngb = 0;
  do
    {
      numngb_inbox = ngb_treefind_variable(&pos[0], h, &startnode);
      for(n = 0; n < numngb_inbox; n++)
	{
	  j = Ngblist[n];
	  dx = pos[0] - P[j].Pos[0];
	  dy = pos[1] - P[j].Pos[1];
	  dz = pos[2] - P[j].Pos[2];
#ifdef PERIODIC			
	  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)
	    {
	      r = sqrt(r2);
	      if(r > 0)
		{
		  u = r * hinv;
		  if(u < 0.5)
		    wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
		  else
		    wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);
		  weight = wk * P[j].Mass / rho;
		  delta_v = sqrt(2 * energy * weight / P[j].Mass);
		  dvx = delta_v * (-dx) / r;
		  dvy = delta_v * (-dy) / r;
		  dvz = delta_v * (-dz) / r;
		  P[j].Vel[0] += dvx;
		  P[j].Vel[1] += dvy;
		  P[j].Vel[2] += dvz;
		  SphP[j].VelPred[0] += dvx;
		  SphP[j].VelPred[1] += dvy;
		  SphP[j].VelPred[2] += dvz;
		}
	    }
	}
    }
  while(startnode >= 0);
}
void kinetic_evaluate ( int target , int mode ) { int j , n ; int startnode , numngb , numngb_inbox ; double h , h2 , hinv , hinv3 ; double rho , weight , wk , energy , delta_v ; double dx , dy , dz , r , r2 , u ; double dvx , dvy , dvz ; MyFloat * pos ; if ( mode == 0 ) { pos = P [ target ] . Pos ; h = PPP [ target ] . Hsml ; rho = SphP [ target ] . Density ; energy = SphP [ target ] . FeedbackEnergy ; } else { pos = KinDataGet [ target ] . Pos ; h = KinDataGet [ target ] . Hsml ; rho = KinDataGet [ target ] . Density ; energy = KinDataGet [ target ] . Energy ; } h2 = h * h ; hinv = 1.0 / h ; #ifndef TWODIMS hinv3 = hinv * hinv * hinv ; #else hinv3 = hinv * hinv / boxSize_Z ; #endif startnode = All . MaxPart ; numngb = 0 ; do { numngb_inbox = ngb_treefind_variable ( & pos [ 0 ] , h , & startnode ) ; for ( n = 0 ; n < numngb_inbox ; n ++ ) { j = Ngblist [ n ] ; dx = pos [ 0 ] - P [ j ] . Pos [ 0 ] ; dy = pos [ 1 ] - P [ j ] . Pos [ 1 ] ; dz = pos [ 2 ] - P [ j ] . Pos [ 2 ] ; #ifdef PERIODIC 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 ) { r = sqrt ( r2 ) ; if ( r > 0 ) { u = r * hinv ; if ( u < 0.5 ) wk = hinv3 * ( KERNEL_COEFF_1 + KERNEL_COEFF_2 * ( u - 1 ) * u * u ) ; else wk = hinv3 * KERNEL_COEFF_5 * ( 1.0 - u ) * ( 1.0 - u ) * ( 1.0 - u ) ; weight = wk * P [ j ] . Mass / rho ; delta_v = sqrt ( 2 * energy * weight / P [ j ] . Mass ) ; dvx = delta_v * ( - dx ) / r ; dvy = delta_v * ( - dy ) / r ; dvz = delta_v * ( - dz ) / r ; P [ j ] . Vel [ 0 ] += dvx ; P [ j ] . Vel [ 1 ] += dvy ; P [ j ] . Vel [ 2 ] += dvz ; SphP [ j ] . VelPred [ 0 ] += dvx ; SphP [ j ] . VelPred [ 1 ] += dvy ; SphP [ j ] . VelPred [ 2 ] += dvz ; } } } } while ( startnode >= 0 ) ; }
This function represents the core of the SPH density computation.
This function represents the core of the SPH density computation .
/*  now find the closest image in the given box size  */