Use floating-point constants and not double-precision ones in the HLLC solver.

src/adiabatic_index.h

@@ -42,6 +42,7 @@
 
 #define hydro_gamma 1.66666666666666667f
 #define hydro_gamma_minus_one 0.66666666666666667f
+#define hydro_gamma_plus_one 2.66666666666666667f
 #define hydro_one_over_gamma_minus_one 1.5f
 #define hydro_gamma_plus_one_over_two_gamma 0.8f
 #define hydro_gamma_minus_one_over_two_gamma 0.2f
@@ -56,6 +57,7 @@
 
 #define hydro_gamma 1.4f
 #define hydro_gamma_minus_one 0.4f
+#define hydro_gamma_plus_one 2.4f
 #define hydro_one_over_gamma_minus_one 2.5f
 #define hydro_gamma_plus_one_over_two_gamma 0.857142857f
 #define hydro_gamma_minus_one_over_two_gamma 0.142857143f
@@ -70,6 +72,7 @@
 
 #define hydro_gamma 1.33333333333333333f
 #define hydro_gamma_minus_one 0.33333333333333333f
+#define hydro_gamma_plus_one 2.33333333333333333f
 #define hydro_one_over_gamma_minus_one 3.f
 #define hydro_gamma_plus_one_over_two_gamma 0.875f
 #define hydro_gamma_minus_one_over_two_gamma 0.125f
@@ -84,6 +87,7 @@
 
 #define hydro_gamma 2.f
 #define hydro_gamma_minus_one 1.f
+#define hydro_gamma_plus_one 3.f
 #define hydro_one_over_gamma_minus_one 1.f
 #define hydro_gamma_plus_one_over_two_gamma 0.75f
 #define hydro_gamma_minus_one_over_two_gamma 0.25f

src/riemann/riemann_exact_isothermal.h

@@ -58,7 +58,7 @@ __attribute__((always_inline)) INLINE static float riemann_fprimeb(float rho,
   if (rho < W[0]) {
     return const_isothermal_soundspeed * W[0] / rho;
   } else {
-    return 0.5 * const_isothermal_soundspeed *
+    return 0.5f * const_isothermal_soundspeed *
            (sqrtf(rho / W[0]) + sqrtf(W[0] / rho)) / rho;
   }
 }

src/riemann/riemann_hllc.h

@@ -40,26 +40,21 @@
 __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
     float *WL, float *WR, float *n, float *vij, float *totflux) {
 
-  float uL, uR, aL, aR;
-  float rhobar, abar, pPVRS, pstar, qL, qR, SL, SR, Sstar;
-  float v2, eL, eR;
-  float UstarL[5], UstarR[5];
-
   /* Handle pure vacuum */
   if (!WL[0] && !WR[0]) {
-    totflux[0] = 0.;
-    totflux[1] = 0.;
-    totflux[2] = 0.;
-    totflux[3] = 0.;
-    totflux[4] = 0.;
+    totflux[0] = 0.f;
+    totflux[1] = 0.f;
+    totflux[2] = 0.f;
+    totflux[3] = 0.f;
+    totflux[4] = 0.f;
     return;
   }
 
   /* STEP 0: obtain velocity in interface frame */
-  uL = WL[1] * n[0] + WL[2] * n[1] + WL[3] * n[2];
-  uR = WR[1] * n[0] + WR[2] * n[1] + WR[3] * n[2];
-  aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
-  aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
+  const float uL = WL[1] * n[0] + WL[2] * n[1] + WL[3] * n[2];
+  const float uR = WR[1] * n[0] + WR[2] * n[1] + WR[3] * n[2];
+  const float aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
+  const float aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
 
   /* Handle vacuum: vacuum does not require iteration and is always exact */
   if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) {
@@ -68,30 +63,30 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
   }
 
   /* STEP 1: pressure estimate */
-  rhobar = 0.5 * (WL[0] + WR[0]);
-  abar = 0.5 * (aL + aR);
-  pPVRS = 0.5 * (WL[4] + WR[4]) - 0.5 * (uR - uL) * rhobar * abar;
-  pstar = max(0., pPVRS);
+  const float rhobar = 0.5f * (WL[0] + WR[0]);
+  const float abar = 0.5f * (aL + aR);
+  const float pPVRS = 0.5f * (WL[4] + WR[4]) - 0.5f * (uR - uL) * rhobar * abar;
+  const float pstar = max(0.f, pPVRS);
 
   /* STEP 2: wave speed estimates
      all these speeds are along the interface normal, since uL and uR are */
-  qL = 1.;
+  float qL = 1.f;
   if (pstar > WL[4]) {
-    qL = sqrtf(1. +
-               0.5 * (hydro_gamma + 1.) / hydro_gamma * (pstar / WL[4] - 1.));
+    qL = sqrtf(1.f +
+               0.5f * hydro_gamma_plus_one * hydro_one_over_gamma * (pstar / WL[4] - 1.f));
   }
-  qR = 1.;
+  float qR = 1.f;
   if (pstar > WR[4]) {
-    qR = sqrtf(1. +
-               0.5 * (hydro_gamma + 1.) / hydro_gamma * (pstar / WR[4] - 1.));
+    qR = sqrtf(1.f +
+               0.5f * hydro_gamma_plus_one * hydro_one_over_gamma * (pstar / WR[4] - 1.f));
   }
-  SL = uL - aL * qL;
-  SR = uR + aR * qR;
-  Sstar = (WR[4] - WL[4] + WL[0] * uL * (SL - uL) - WR[0] * uR * (SR - uR)) /
-          (WL[0] * (SL - uL) - WR[0] * (SR - uR));
+  const float SL = uL - aL * qL;
+  const float SR = uR + aR * qR;
+  const float Sstar = (WR[4] - WL[4] + WL[0] * uL * (SL - uL) - WR[0] * uR * (SR - uR)) /
+    (WL[0] * (SL - uL) - WR[0] * (SR - uR));
 
   /* STEP 3: HLLC flux in a frame moving with the interface velocity */
-  if (Sstar >= 0.) {
+  if (Sstar >= 0.f) {
     /* flux FL */
     totflux[0] = WL[0] * uL;
     /* these are the actual correct fluxes in the boosted lab frame
@@ -99,15 +94,18 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
     totflux[1] = WL[0] * WL[1] * uL + WL[4] * n[0];
     totflux[2] = WL[0] * WL[2] * uL + WL[4] * n[1];
     totflux[3] = WL[0] * WL[3] * uL + WL[4] * n[2];
-    v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3];
-    eL = WL[4] / hydro_gamma_minus_one / WL[0] + 0.5 * v2;
+    const float v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3];
+    const float eL = WL[4] * hydro_one_over_gamma_minus_one / WL[0] + 0.5f * v2;
     totflux[4] = WL[0] * eL * uL + WL[4] * uL;
-    if (SL < 0.) {
+    if (SL < 0.f) {
+
+      float UstarL[5];
+
       /* add flux FstarL */
-      UstarL[0] = 1.;
+      UstarL[0] = 1.f;
       /* we need UstarL in the lab frame:
-         subtract the velocity in the interface frame from the lab frame
-         velocity and then add Sstar in interface frame */
+       * subtract the velocity in the interface frame from the lab frame
+       * velocity and then add Sstar in interface frame */
       UstarL[1] = WL[1] + (Sstar - uL) * n[0];
       UstarL[2] = WL[2] + (Sstar - uL) * n[1];
       UstarL[3] = WL[3] + (Sstar - uL) * n[2];
@@ -129,12 +127,19 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
     totflux[1] = WR[0] * WR[1] * uR + WR[4] * n[0];
     totflux[2] = WR[0] * WR[2] * uR + WR[4] * n[1];
     totflux[3] = WR[0] * WR[3] * uR + WR[4] * n[2];
-    v2 = WR[1] * WR[1] + WR[2] * WR[2] + WR[3] * WR[3];
-    eR = WR[4] / hydro_gamma_minus_one / WR[0] + 0.5 * v2;
+    const float v2 = WR[1] * WR[1] + WR[2] * WR[2] + WR[3] * WR[3];
+    const float eR = WR[4] * hydro_one_over_gamma_minus_one / WR[0] + 0.5f * v2;
     totflux[4] = WR[0] * eR * uR + WR[4] * uR;
-    if (SR > 0.) {
+    if (SR > 0.f) {
+
+      float UstarR[5];
+
       /* add flux FstarR */
-      UstarR[0] = 1.;
+      UstarR[0] = 1.f;
+
+      /* we need UstarR in the lab frame:
+       * subtract the velocity in the interface frame from the lab frame
+       * velocity and then add Sstar in interface frame */
       UstarR[1] = WR[1] + (Sstar - uR) * n[0];
       UstarR[2] = WR[2] + (Sstar - uR) * n[1];
       UstarR[3] = WR[3] + (Sstar - uR) * n[2];
@@ -152,16 +157,16 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
     }
   }
 
-  /* deboost to lab frame
-     we add the flux contribution due to the movement of the interface
-     the density flux is unchanged
+  /* deboost to lab frame we add the flux contribution due to the 
+     movement of the interface the density flux is unchanged
      we add the extra velocity flux due to the absolute motion of the fluid
      similarly, we need to add the energy fluxes due to the absolute motion */
-  v2 = vij[0] * vij[0] + vij[1] * vij[1] + vij[2] * vij[2];
-  // order is important: we first use the momentum fluxes to update the energy
-  // flux and then de-boost the momentum fluxes!
+  const float v2 = vij[0] * vij[0] + vij[1] * vij[1] + vij[2] * vij[2];
+
+  /* order is important: we first use the momentum fluxes to update the energy
+     flux and then de-boost the momentum fluxes! */
   totflux[4] += vij[0] * totflux[1] + vij[1] * totflux[2] +
-                vij[2] * totflux[3] + 0.5 * v2 * totflux[0];
+                vij[2] * totflux[3] + 0.5f * v2 * totflux[0];
   totflux[1] += vij[0] * totflux[0];
   totflux[2] += vij[1] * totflux[0];
   totflux[3] += vij[2] * totflux[0];