Merge branch 'SmoothLengthFix' into 'master'

Removed an extra kernel_igamma from the rho_dh calculation

There was an extra kernel_igamma in the calculation of rho_sh which was causing an offset to equivalent Gadget-2 runs. 

See merge request !162

src/hydro/Default/hydro.h

@@ -93,7 +93,7 @@ __attribute__((always_inline))
 
   /* Final operation on the density (add self-contribution). */
   p->rho += p->mass * kernel_root;
-  p->rho_dh -= 3.0f * p->mass * kernel_root * kernel_igamma;
+  p->rho_dh -= 3.0f * p->mass * kernel_root;
   p->density.wcount += kernel_root;
 
   /* Finish the calculation by inserting the missing h-factors */
@@ -101,6 +101,11 @@ __attribute__((always_inline))
   p->rho_dh *= ih4;
   p->density.wcount *= (4.0f / 3.0f * M_PI * kernel_gamma3);
   p->density.wcount_dh *= ih * (4.0f / 3.0f * M_PI * kernel_gamma4);
+
+  const float irho = 1.f / p->rho;
+
+  /* Compute the derivative term */
+  p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
 }
 
 /**

src/hydro/Gadget2/hydro.h

@@ -94,7 +94,7 @@ __attribute__((always_inline))
 
   /* Final operation on the density (add self-contribution). */
   p->rho += p->mass * kernel_root;
-  p->rho_dh -= 3.0f * p->mass * kernel_root * kernel_igamma;
+  p->rho_dh -= 3.0f * p->mass * kernel_root;
   p->density.wcount += kernel_root;
 
   /* Finish the calculation by inserting the missing h-factors */

src/hydro/Gadget2/hydro_iact.h

@@ -65,7 +65,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   /* Compute contribution to the density */
   pi->rho += mj * wi;
-  pi->rho_dh -= mj * kernel_igamma * (3.f * wi + ui * wi_dx);
+  pi->rho_dh -= mj * (3.f * wi + ui * wi_dx);
 
   /* Compute contribution to the number of neighbours */
   pi->density.wcount += wi;
@@ -78,7 +78,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   /* Compute contribution to the density */
   pj->rho += mi * wj;
-  pj->rho_dh -= mi * kernel_igamma * (3.f * wj + uj * wj_dx);
+  pj->rho_dh -= mi * (3.f * wj + uj * wj_dx);
 
   /* Compute contribution to the number of neighbours */
   pj->density.wcount += wj;
@@ -134,7 +134,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
 
   /* Compute contribution to the density */
   pi->rho += mj * wi;
-  pi->rho_dh -= mj * kernel_igamma * (3.f * wi + u * wi_dx);
+  pi->rho_dh -= mj * (3.f * wi + u * wi_dx);
 
   /* Compute contribution to the number of neighbours */
   pi->density.wcount += wi;

src/hydro/Minimal/hydro.h

@@ -94,7 +94,7 @@ __attribute__((always_inline))
 
   /* Final operation on the density (add self-contribution). */
   p->rho += p->mass * kernel_root;
-  p->rho_dh -= 3.0f * p->mass * kernel_root * kernel_igamma;
+  p->rho_dh -= 3.0f * p->mass * kernel_root;
   p->density.wcount += kernel_root;
 
   /* Finish the calculation by inserting the missing h-factors */