Accelerations and time-steps now match Gadget-2

src/const.h

@@ -38,7 +38,7 @@
   1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
 
 /* Time integration constants. */
-#define const_cfl 0.3f
+#define const_cfl 0.2f
 #define const_ln_max_h_change                                           \
   0.231111721f /* Particle can't change volume by more than a factor of \
                   2=1.26^3 over one time step */
@@ -47,7 +47,7 @@
 /* Neighbour search constants. */
 #define const_eta_kernel \
   1.2349f /* Corresponds to 48 ngbs with the cubic spline kernel */
-#define const_delta_nwneigh 1.f
+#define const_delta_nwneigh 0.1f
 #define const_smoothing_max_iter 30
 #define CUBIC_SPLINE_KERNEL
 

src/debug.c

@@ -50,7 +50,7 @@ void printParticle(struct part *parts, long long int id, int N) {
           "v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e],\n h=%.3e, "
           "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, dS/dt=%.3e,\n"
 	  "divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e]  \n "
-	  "t_begin=%.3e, t_end=%.3e\n",
+	  "v_sig=%e t_begin=%.3e, t_end=%.3e\n",
           i, parts[i].id, parts[i].x[0], parts[i].x[1], parts[i].x[2],
           parts[i].v[0], parts[i].v[1], parts[i].v[2], parts[i].a[0],
           parts[i].a[1], parts[i].a[2], 2.*parts[i].h,
@@ -65,6 +65,7 @@ void printParticle(struct part *parts, long long int id, int N) {
 	  parts[i].rot_v[0],
 	  parts[i].rot_v[1],
 	  parts[i].rot_v[2],
+	  parts[i].v_sig,
           parts[i].t_begin, parts[i].t_end);
       found = 1;
     }

src/hydro/Gadget2/hydro.h

@@ -27,21 +27,16 @@
 __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
     struct part* p, struct xpart* xp) {
 
-  /* CFL condition */
-  const float dt_cfl = const_cfl * p->h / p->v_sig;
-
-  /* /\* Limit change in h *\/ */
-  /* float dt_h_change = (p->h_dt != 0.0f) */
-  /*                         ? fabsf(const_ln_max_h_change * p->h / p->h_dt) */
-  /*                         : FLT_MAX; */
+  /* Acceleration */
+  float ac = sqrtf(p->a[0]*p->a[0] + p->a[1]*p->a[1] + p->a[2]*p->a[2]);
+  ac = fmaxf(ac, 1e-30);
 
-  /* /\* Limit change in u *\/ */
-  /* float dt_u_change = (p->force.u_dt != 0.0f) */
-  /*                         ? fabsf(const_max_u_change * p->u / p->force.u_dt) */
-  /*                         : FLT_MAX; */
+  const float dt_accel = sqrtf(2.f);
+  
+  /* CFL condition */
+  const float dt_cfl = 2.f * const_cfl * p->h / p->v_sig;
 
-  //  return fminf(dt_cfl, fminf(dt_h_change, dt_u_change));
-  return dt_cfl; 
+  return fminf(dt_cfl, dt_accel); 
 }
 
 
@@ -176,8 +171,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(struc
   /* Reset the time derivatives. */
   p->entropy_dt = 0.0f;
   
-  /* Reset minimal signal velocity */
-  p->v_sig = FLT_MAX;
+  /* Reset maximal signal velocity */
+  p->v_sig = 0.0f;
 }
 
 
@@ -228,6 +223,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(struct part* p
  */
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(struct part* p) {
 
-  p->entropy = (const_hydro_gamma - 1.f) * p->entropy * pow(p->rho, -(const_hydro_gamma - 1.f));
+  p->entropy = (const_hydro_gamma - 1.f) * p->entropy * powf(p->rho, -(const_hydro_gamma - 1.f));
 
 }

src/hydro/Gadget2/hydro_iact.h

@@ -113,13 +113,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   pj->rot_v[0] += facj * curlvr[0];
   pj->rot_v[1] += facj * curlvr[1];
   pj->rot_v[2] += facj * curlvr[2];
-
-  /* if(pi->id == 1000) */
-  /*   message("Interacting with %lld. r=%f\n", pj->id, r); */
-
-  /* if(pj->id == 1000) */
-  /*   message("Interacting with %lld. r=%f\n", pi->id, r); */
-
 }
 
 
@@ -139,8 +132,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   /* Get r and r inverse. */
   const float r = sqrtf(r2);
   const float ri = 1.0f / r;
-
-  if(pi->id == 1000 && pj->id == 1103) hi = 0.2234976 / 2.f;
   
   /* Compute the kernel function */
   const float h_inv = 1.0f / hi;
@@ -167,7 +158,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
   pi->div_v -= fac * dvdr;
 
-  if(pi->id == 1000 && pj->id == 1103)
+  if(pi->id == 1000 && pj->id == 1203)
     message("Interacting with %lld. r=%e hi=%e u=%e W=%e dW/dx=%e dh_drho1=%e dh_drho2=%e\n fac=%e dvdr=%e",
 	    pj->id,
 	    r,
@@ -260,6 +251,23 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   /* Eventually got the acceleration */
   const float acc = visc_term + sph_term;
 
+
+  if(pi->id == 1000 && pj->id == 1203)
+    message("Interacting with %lld. r=%e hi=%e hj=%e dWi/dx=%e dWj/dx=%3e dvdr=%e visc=%e sph=%e",
+	    pj->id,
+	    r,
+	    2*hi,
+	    2*hj,
+	    wi_dr,
+	    wj_dr,
+	    dvdr,
+	    visc_term,
+	    sph_term
+	    );
+  if(pi->id == 1203 && pj->id == 1000)
+    message("oO");
+
+  
   /* Use the force Luke ! */
   pi->a[0] -= acc * dx[0];
   pi->a[1] -= acc * dx[1];
@@ -274,8 +282,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   pj->v_sig = fmaxf(pj->v_sig, v_sig) ;
   
   /* Change in entropy */
-  pi->entropy_dt += 0.5f * visc_term * dvdr;
-  pj->entropy_dt -= 0.5f * visc_term * dvdr;
+  pi->entropy_dt -= 0.5f * visc_term * dvdr;
+  pj->entropy_dt += 0.5f * visc_term * dvdr;
 }
 
 
@@ -345,7 +353,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Eventually got the acceleration */
   const float acc = visc_term + sph_term;
-
+  
   /* Use the force Luke ! */
   pi->a[0] -= acc * dx[0];
   pi->a[1] -= acc * dx[1];
@@ -355,7 +363,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   pi->v_sig = fmaxf(pi->v_sig, v_sig) ;
   
   /* Change in entropy */
-  pi->entropy_dt += 0.5f * visc_term * dvdr;
+  pi->entropy_dt -= 0.5f * visc_term * dvdr;
 }
 
 

src/runner.c

@@ -725,7 +725,7 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
       p->x[0] += xp->v_full[0] * dt;
       p->x[1] += xp->v_full[1] * dt;
       p->x[2] += xp->v_full[2] * dt;
-
+      
       /* Predict velocities */
       p->v[0] += p->a[0] * dt;
       p->v[1] += p->a[1] * dt;
@@ -860,6 +860,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
 	  
 	  new_dt = fminf(new_dt_hydro, new_dt_grav);
 
+	  if(p->id == 1000)
+	    message("1000   dt_hydro=%e", new_dt_hydro);
+
+	  if(p->id == 515050)
+	    message("515050 dt_hydro=%e", new_dt_hydro);
+
+	  
 	  /* Recover the current timestep */
 	  const float current_dt = p->t_end - p->t_begin;