Re-ordered some operations in the scalar versions of SPH schemes to take...

Re-ordered some operations in the scalar versions of SPH schemes to take advantage of the faster rsqrt calls.

src/hydro/Gadget2/hydro_iact.h

@@ -59,9 +59,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   const float mi = pi->mass;
   const float mj = pj->mass;
 
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Compute the kernel function for pi */
   const float hi_inv = 1.f / hi;
@@ -148,9 +148,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   /* Get the masses. */
   const float mj = pj->mass;
 
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Compute the kernel function */
   const float hi_inv = 1.0f / hi;
@@ -440,8 +440,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Get some values in local variables. */
   const float mi = pi->mass;
@@ -559,8 +560,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Get some values in local variables. */
   // const float mi = pi->mass;

src/hydro/Minimal/hydro_iact.h

@@ -53,7 +53,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   float wi, wj, wi_dx, wj_dx;
 
-  const float r = sqrtf(r2);
+  /* Get r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Get the masses. */
   const float mi = pi->mass;
@@ -101,8 +103,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   /* Get the masses. */
   const float mj = pj->mass;
 
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
+  /* Get r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   const float h_inv = 1.f / hi;
   const float ui = r * h_inv;
@@ -134,8 +137,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and r inverse. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Recover some data */
   const float mi = pi->mass;
@@ -246,8 +250,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and r inverse. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Recover some data */
   // const float mi = pi->mass;

src/hydro/PressureEntropy/hydro_iact.h

@@ -54,9 +54,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   const float mi = pi->mass;
   const float mj = pj->mass;
 
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Compute the kernel function for pi */
   const float hi_inv = 1.f / hi;
@@ -142,9 +142,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   /* Get the masses. */
   const float mj = pj->mass;
 
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
-  const float ri = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Compute the kernel function */
   const float h_inv = 1.0f / hi;
@@ -164,7 +164,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   pi->density.pressure_dh -=
       mj * pj->entropy_one_over_gamma * (hydro_dimension * wi + ui * wi_dx);
 
-  const float fac = mj * wi_dx * ri;
+  const float fac = mj * wi_dx * r_inv;
 
   /* Compute dv dot r */
   dv[0] = pi->v[0] - pj->v[0];
@@ -205,8 +205,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r . */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Get some values in local variables. */
   const float mi = pi->mass;
@@ -318,8 +319,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Get some values in local variables. */
   // const float mi = pi->mass;