Velocity curl and divergence now correct as well.

src/debug.c

@@ -48,8 +48,8 @@ void printParticle(struct part *parts, long long int id, int N) {
       printf(
           "## Particle[%d]:\n id=%lld, x=[%.3e,%.3e,%.3e], "
           "v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e],\n h=%.3e, "
-          "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%3.e, S=%.3e,\n divV=%.3e, "
-	  " curlV=%.3e rotV=[%.3e,%.3e,%.3e]  \n "
+          "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",
           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],
@@ -59,6 +59,7 @@ void printParticle(struct part *parts, long long int id, int N) {
 	  parts[i].rho,
 	  parts[i].pressure,
 	  parts[i].entropy,
+	  parts[i].entropy_dt,
 	  parts[i].div_v,
 	  parts[i].curl_v,
 	  parts[i].rot_v[0],

src/hydro/Gadget2/hydro.h

@@ -98,16 +98,21 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(struct part*
   p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh / p->rho);
 
   /* Finish calculation of the velocity curl */
+  p->rot_v[0] *= ih4;
+  p->rot_v[1] *= ih4;
+  p->rot_v[2] *= ih4;
+
+  /* And compute the norm of the curl */
   p->curl_v = sqrtf(p->rot_v[0] * p->rot_v[0] +
 		    p->rot_v[1] * p->rot_v[1] +
 		    p->rot_v[2] * p->rot_v[2]) / p->rho;
   
   /* Finish calculation of the velocity divergence */
-  p->div_v /= p->rho;
+  p->div_v *= ih4 / p->rho;
   
   /* Compute the pressure */
   const float dt = time - 0.5f*(p->t_begin + p->t_end);
-  p->pressure = (p->entropy + p->entropy_dt * dt) * pow(p->rho, const_hydro_gamma);
+  p->pressure = (p->entropy + p->entropy_dt * dt) * powf(p->rho, const_hydro_gamma);
 
 }
 

src/hydro/Gadget2/hydro_iact.h

@@ -98,8 +98,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   dv[2] = pi->v[2] - pj->v[2];
   const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
 
-  pi->div_v -= faci * dvdr;
-  pj->div_v -= facj * dvdr;
+  pi->div_v += faci * dvdr;
+  pj->div_v += facj * dvdr;
 
   /* Compute dv cross r */
   curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
@@ -168,7 +168,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   pi->div_v -= fac * dvdr;
 
   if(pi->id == 1000 && pj->id == 1103)
-    message("Interacting with %lld. r=%e hi=%e u=%e W=%e dW/dx=%e dh_drho1=%e dh_drho2=%e %e\n fac=%e dvdr=%e",
+    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,
 	    hi,
@@ -177,7 +177,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
 	    wi_dx * ih4,
 	    -mj * (3.f * kernel_igamma * wi) * ih4,
 	    -mj * u * wi_dx * kernel_igamma * ih4,
-	    fac,
+	    fac * ih4,
 	    dvdr
 	    );