Merge branch 'gizmo_mfm_clean' into 'master'

Gizmo mfm clean

See merge request !588

examples/SineWavePotential_1D/makeIC.py

@@ -31,12 +31,12 @@ cs2 = 2.*uconst/3.
 A = 10.
 
 fileName = "sineWavePotential.hdf5"
-numPart = 100
+numPart = 1000
 boxSize = 1.
 
 coords = np.zeros((numPart, 3))
 v = np.zeros((numPart, 3))
-m = np.zeros(numPart) + 1.
+m = np.zeros(numPart) + 1000. / numPart
 h = np.zeros(numPart) + 2./numPart
 u = np.zeros(numPart) + uconst
 ids = np.arange(numPart, dtype = 'L')

examples/SineWavePotential_1D/plotSolution.py

@@ -23,8 +23,13 @@
 import numpy as np
 import h5py
 import sys
+import matplotlib
+matplotlib.use("Agg")
 import pylab as pl
 
+pl.rcParams["figure.figsize"] = (12, 10)
+pl.rcParams["text.usetex"] = True
+
 # these should be the same as in makeIC.py
 uconst = 20.2615290634
 cs2 = 2.*uconst/3.
@@ -39,15 +44,20 @@ fileName = sys.argv[1]
 file = h5py.File(fileName, 'r')
 coords = np.array(file["/PartType0/Coordinates"])
 rho = np.array(file["/PartType0/Density"])
+P = np.array(file["/PartType0/Pressure"])
 u = np.array(file["/PartType0/InternalEnergy"])
-agrav = np.array(file["/PartType0/GravAcceleration"])
 m = np.array(file["/PartType0/Masses"])
+vs = np.array(file["/PartType0/Velocities"])
 ids = np.array(file["/PartType0/ParticleIDs"])
 
 x = np.linspace(0., 1., 1000)
 rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
 
-P = cs2*rho
+a = A * np.sin(2. * np.pi * x)
+
+apart = A * np.sin(2. * np.pi * coords[:,0])
+tkin = -0.5 * np.dot(apart, coords[:,0])
+print tkin, 0.5 * np.dot(m, vs[:,0]**2)
 
 ids_reverse = np.argsort(ids)
 
@@ -65,13 +75,38 @@ for i in range(len(P)):
     corr = 1.
   idxp1 = ids_reverse[ip1]
   idxm1 = ids_reverse[im1]
-  gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0])
+  gradP[i] = (P[idxp1] - P[idxm1])/(coords[idxp1,0] - coords[idxm1,0])
+
+fig, ax = pl.subplots(2, 2, sharex = True)
+
+ax[0][0].plot(coords[:,0], rho, ".", label = "gizmo-mfm")
+ax[0][0].plot(x, rho_x, "-", label = "stable solution")
+ax[0][0].set_ylabel("$\\rho{}$ (code units)")
+ax[0][0].legend(loc = "best")
+
+ax[0][1].plot(x, a, "-", label = "$\\nabla{}\\Phi{}$ external")
+ax[0][1].plot(coords[:,0], gradP/rho, ".",
+              label = "$\\nabla{}P/\\rho{}$ gizmo-mfm")
+ax[0][1].set_ylabel("force (code units)")
+ax[0][1].legend(loc = "best")
+
+ax[1][0].axhline(y = uconst, label = "isothermal $u$", color = "k",
+                 linestyle = "--")
+ax[1][0].plot(coords[:,0], u, ".", label = "gizmo-mfm")
+ax[1][0].set_ylabel("$u$ (code units)")
+ax[1][0].set_xlabel("$x$ (code units)")
+ax[1][0].legend(loc = "best")
+
+#ax[1][1].plot(coords[:,0], m, "y.")
+#ax[1][1].set_ylabel("$m_i$ (code units)")
+#ax[1][1].set_xlabel("$x$ (code units)")
 
-fig, ax = pl.subplots(2, 2)
+ax[1][1].axhline(y = 0., label = "target", color = "k",
+                 linestyle = "--")
+ax[1][1].plot(coords[:,0], vs[:,0], ".", label = "gizmo-mfm")
+ax[1][1].set_ylabel("$v_x$ (code units)")
+ax[1][1].set_xlabel("$x$ (code units)")
+ax[1][1].legend(loc = "best")
 
-ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
-ax[0][0].plot(x, rho_x, "g-")
-ax[0][1].plot(coords[:,0], gradP/rho, "b.")
-ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5)
-ax[1][1].plot(coords[:,0], m, "y.")
+pl.tight_layout()
 pl.savefig("{fileName}.png".format(fileName = fileName[:-5]))

examples/SineWavePotential_1D/sineWavePotential.yml

@@ -10,7 +10,7 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   10.   # The end time of the simulation (in internal units).
-  dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
+  dt_min:     1e-8  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the conserved quantities statistics
@@ -36,3 +36,6 @@ InitialConditions:
 SineWavePotential:
   amplitude: 10.
   timestep_limit: 1.
+
+EoS:
+  isothermal_internal_energy: 20.2615290634

src/Makefile.am

@@ -106,7 +106,6 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  hydro/GizmoMFM/hydro_slope_limiters_face.h \
                  hydro/GizmoMFM/hydro_slope_limiters.h \
                  hydro/GizmoMFM/hydro_unphysical.h \
-                 hydro/GizmoMFM/hydro_velocities.h \
                  hydro/Shadowswift/hydro_debug.h \
                  hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
                  hydro/Shadowswift/hydro_iact.h \

src/hydro/GizmoMFM/hydro_debug.h

@@ -27,8 +27,6 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "a=[%.3e,%.3e,%.3e], "
       "h=%.3e, "
       "time_bin=%d, "
-      "primitives={"
-      "v=[%.3e,%.3e,%.3e], "
       "rho=%.3e, "
       "P=%.3e, "
       "gradients={"
@@ -39,7 +37,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "rho=[%.3e,%.3e], "
       "v=[[%.3e,%.3e],[%.3e,%.3e],[%.3e,%.3e]], "
       "P=[%.3e,%.3e], "
-      "maxr=%.3e}}, "
+      "maxr=%.3e}, "
       "conserved={"
       "momentum=[%.3e,%.3e,%.3e], "
       "mass=%.3e, "
@@ -53,24 +51,18 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "wcount_dh=%.3e, "
       "wcount=%.3e}\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
-      p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->primitives.v[0],
-      p->primitives.v[1], p->primitives.v[2], p->primitives.rho,
-      p->primitives.P, p->primitives.gradients.rho[0],
-      p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
-      p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
-      p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
-      p->primitives.gradients.v[1][1], p->primitives.gradients.v[1][2],
-      p->primitives.gradients.v[2][0], p->primitives.gradients.v[2][1],
-      p->primitives.gradients.v[2][2], p->primitives.gradients.P[0],
-      p->primitives.gradients.P[1], p->primitives.gradients.P[2],
-      p->primitives.limiter.rho[0], p->primitives.limiter.rho[1],
-      p->primitives.limiter.v[0][0], p->primitives.limiter.v[0][1],
-      p->primitives.limiter.v[1][0], p->primitives.limiter.v[1][1],
-      p->primitives.limiter.v[2][0], p->primitives.limiter.v[2][1],
-      p->primitives.limiter.P[0], p->primitives.limiter.P[1],
-      p->primitives.limiter.maxr, p->conserved.momentum[0],
-      p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.mass,
-      p->conserved.energy, p->geometry.volume, p->geometry.matrix_E[0][0],
+      p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->rho, p->P,
+      p->gradients.rho[0], p->gradients.rho[1], p->gradients.rho[2],
+      p->gradients.v[0][0], p->gradients.v[0][1], p->gradients.v[0][2],
+      p->gradients.v[1][0], p->gradients.v[1][1], p->gradients.v[1][2],
+      p->gradients.v[2][0], p->gradients.v[2][1], p->gradients.v[2][2],
+      p->gradients.P[0], p->gradients.P[1], p->gradients.P[2],
+      p->limiter.rho[0], p->limiter.rho[1], p->limiter.v[0][0],
+      p->limiter.v[0][1], p->limiter.v[1][0], p->limiter.v[1][1],
+      p->limiter.v[2][0], p->limiter.v[2][1], p->limiter.P[0], p->limiter.P[1],
+      p->limiter.maxr, p->conserved.momentum[0], p->conserved.momentum[1],
+      p->conserved.momentum[2], p->conserved.mass, p->conserved.energy,
+      p->geometry.volume, p->geometry.matrix_E[0][0],
       p->geometry.matrix_E[0][1], p->geometry.matrix_E[0][2],
       p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
       p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],

src/hydro/GizmoMFM/hydro_gradients.h

@@ -101,38 +101,28 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
   const float xij_j[3] = {xij_i[0] + dx[0], xij_i[1] + dx[1], xij_i[2] + dx[2]};
 
   float dWi[5];
-  dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +
-           pi->primitives.gradients.rho[1] * xij_i[1] +
-           pi->primitives.gradients.rho[2] * xij_i[2];
-  dWi[1] = pi->primitives.gradients.v[0][0] * xij_i[0] +
-           pi->primitives.gradients.v[0][1] * xij_i[1] +
-           pi->primitives.gradients.v[0][2] * xij_i[2];
-  dWi[2] = pi->primitives.gradients.v[1][0] * xij_i[0] +
-           pi->primitives.gradients.v[1][1] * xij_i[1] +
-           pi->primitives.gradients.v[1][2] * xij_i[2];
-  dWi[3] = pi->primitives.gradients.v[2][0] * xij_i[0] +
-           pi->primitives.gradients.v[2][1] * xij_i[1] +
-           pi->primitives.gradients.v[2][2] * xij_i[2];
-  dWi[4] = pi->primitives.gradients.P[0] * xij_i[0] +
-           pi->primitives.gradients.P[1] * xij_i[1] +
-           pi->primitives.gradients.P[2] * xij_i[2];
+  dWi[0] = pi->gradients.rho[0] * xij_i[0] + pi->gradients.rho[1] * xij_i[1] +
+           pi->gradients.rho[2] * xij_i[2];
+  dWi[1] = pi->gradients.v[0][0] * xij_i[0] + pi->gradients.v[0][1] * xij_i[1] +
+           pi->gradients.v[0][2] * xij_i[2];
+  dWi[2] = pi->gradients.v[1][0] * xij_i[0] + pi->gradients.v[1][1] * xij_i[1] +
+           pi->gradients.v[1][2] * xij_i[2];
+  dWi[3] = pi->gradients.v[2][0] * xij_i[0] + pi->gradients.v[2][1] * xij_i[1] +
+           pi->gradients.v[2][2] * xij_i[2];
+  dWi[4] = pi->gradients.P[0] * xij_i[0] + pi->gradients.P[1] * xij_i[1] +
+           pi->gradients.P[2] * xij_i[2];
 
   float dWj[5];
-  dWj[0] = pj->primitives.gradients.rho[0] * xij_j[0] +
-           pj->primitives.gradients.rho[1] * xij_j[1] +
-           pj->primitives.gradients.rho[2] * xij_j[2];
-  dWj[1] = pj->primitives.gradients.v[0][0] * xij_j[0] +
-           pj->primitives.gradients.v[0][1] * xij_j[1] +
-           pj->primitives.gradients.v[0][2] * xij_j[2];
-  dWj[2] = pj->primitives.gradients.v[1][0] * xij_j[0] +
-           pj->primitives.gradients.v[1][1] * xij_j[1] +
-           pj->primitives.gradients.v[1][2] * xij_j[2];
-  dWj[3] = pj->primitives.gradients.v[2][0] * xij_j[0] +
-           pj->primitives.gradients.v[2][1] * xij_j[1] +
-           pj->primitives.gradients.v[2][2] * xij_j[2];
-  dWj[4] = pj->primitives.gradients.P[0] * xij_j[0] +
-           pj->primitives.gradients.P[1] * xij_j[1] +
-           pj->primitives.gradients.P[2] * xij_j[2];
+  dWj[0] = pj->gradients.rho[0] * xij_j[0] + pj->gradients.rho[1] * xij_j[1] +
+           pj->gradients.rho[2] * xij_j[2];
+  dWj[1] = pj->gradients.v[0][0] * xij_j[0] + pj->gradients.v[0][1] * xij_j[1] +
+           pj->gradients.v[0][2] * xij_j[2];
+  dWj[2] = pj->gradients.v[1][0] * xij_j[0] + pj->gradients.v[1][1] * xij_j[1] +
+           pj->gradients.v[1][2] * xij_j[2];
+  dWj[3] = pj->gradients.v[2][0] * xij_j[0] + pj->gradients.v[2][1] * xij_j[1] +
+           pj->gradients.v[2][2] * xij_j[2];
+  dWj[4] = pj->gradients.P[0] * xij_j[0] + pj->gradients.P[1] * xij_j[1] +
+           pj->gradients.P[2] * xij_j[2];
 
   /* Apply the slope limiter at this interface */
   hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);

src/hydro/GizmoMFM/hydro_gradients_sph.h

@@ -28,24 +28,24 @@
 __attribute__((always_inline)) INLINE static void hydro_gradients_init(
     struct part *p) {
 
-  p->primitives.gradients.rho[0] = 0.0f;
-  p->primitives.gradients.rho[1] = 0.0f;
-  p->primitives.gradients.rho[2] = 0.0f;
+  p->gradients.rho[0] = 0.0f;
+  p->gradients.rho[1] = 0.0f;
+  p->gradients.rho[2] = 0.0f;
 
-  p->primitives.gradients.v[0][0] = 0.0f;
-  p->primitives.gradients.v[0][1] = 0.0f;
-  p->primitives.gradients.v[0][2] = 0.0f;
+  p->gradients.v[0][0] = 0.0f;
+  p->gradients.v[0][1] = 0.0f;
+  p->gradients.v[0][2] = 0.0f;
 
-  p->primitives.gradients.v[1][0] = 0.0f;
-  p->primitives.gradients.v[1][1] = 0.0f;
-  p->primitives.gradients.v[1][2] = 0.0f;
-  p->primitives.gradients.v[2][0] = 0.0f;
-  p->primitives.gradients.v[2][1] = 0.0f;
-  p->primitives.gradients.v[2][2] = 0.0f;
+  p->gradients.v[1][0] = 0.0f;
+  p->gradients.v[1][1] = 0.0f;
+  p->gradients.v[1][2] = 0.0f;
+  p->gradients.v[2][0] = 0.0f;
+  p->gradients.v[2][1] = 0.0f;
+  p->gradients.v[2][2] = 0.0f;
 
-  p->primitives.gradients.P[0] = 0.0f;
-  p->primitives.gradients.P[1] = 0.0f;
-  p->primitives.gradients.P[2] = 0.0f;
+  p->gradients.P[0] = 0.0f;
+  p->gradients.P[1] = 0.0f;
+  p->gradients.P[2] = 0.0f;
 
   hydro_slope_limit_cell_init(p);
 }
@@ -64,7 +64,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
     float r2, const float *dx, float hi, float hj, struct part *restrict pi,
     struct part *restrict pj) {
 
-  const float r_inv = 1.f / sqrtf(r2);
+  const float r_inv = 1.0f / sqrtf(r2);
   const float r = r2 * r_inv;
 
   float wi, wi_dx;
@@ -72,41 +72,29 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
   const float xi = r * hi_inv;
   kernel_deval(xi, &wi, &wi_dx);
 
-  /* very basic gradient estimate */
-  pi->primitives.gradients.rho[0] -=
-      wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-  pi->primitives.gradients.rho[1] -=
-      wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-  pi->primitives.gradients.rho[2] -=
-      wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-
-  pi->primitives.gradients.v[0][0] -=
-      wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-  pi->primitives.gradients.v[0][1] -=
-      wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-  pi->primitives.gradients.v[0][2] -=
-      wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-
-  pi->primitives.gradients.v[1][0] -=
-      wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-  pi->primitives.gradients.v[1][1] -=
-      wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-  pi->primitives.gradients.v[1][2] -=
-      wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-
-  pi->primitives.gradients.v[2][0] -=
-      wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-  pi->primitives.gradients.v[2][1] -=
-      wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-  pi->primitives.gradients.v[2][2] -=
-      wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-
-  pi->primitives.gradients.P[0] -=
-      wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
-  pi->primitives.gradients.P[1] -=
-      wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
-  pi->primitives.gradients.P[2] -=
-      wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+  const float dW[5] = {pi->rho - pj->rho, pi->v[0] - pj->v[0],
+                       pi->v[1] - pj->v[1], pi->v[2] - pj->v[2], pi->P - pj->P};
+
+  const float normi = wi_dx * r_inv;
+  const float nidx[3] = {normi * dx[0], normi * dx[1], normi * dx[2]};
+
+  pi->gradients.rho[0] -= dW[0] * nidx[0];
+  pi->gradients.rho[1] -= dW[0] * nidx[1];
+  pi->gradients.rho[2] -= dW[0] * nidx[2];
+
+  pi->gradients.v[0][0] -= dW[1] * nidx[0];
+  pi->gradients.v[0][1] -= dW[1] * nidx[1];
+  pi->gradients.v[0][2] -= dW[1] * nidx[2];
+  pi->gradients.v[1][0] -= dW[2] * nidx[0];
+  pi->gradients.v[1][1] -= dW[2] * nidx[1];
+  pi->gradients.v[1][2] -= dW[2] * nidx[2];
+  pi->gradients.v[2][0] -= dW[3] * nidx[0];
+  pi->gradients.v[2][1] -= dW[3] * nidx[1];
+  pi->gradients.v[2][2] -= dW[3] * nidx[2];
+
+  pi->gradients.P[0] -= dW[4] * nidx[0];
+  pi->gradients.P[1] -= dW[4] * nidx[1];
+  pi->gradients.P[2] -= dW[4] * nidx[2];
 
   hydro_slope_limit_cell_collect(pi, pj, r);
 
@@ -115,40 +103,27 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
   const float xj = r * hj_inv;
   kernel_deval(xj, &wj, &wj_dx);
 
+  const float normj = wj_dx * r_inv;
+  const float njdx[3] = {normj * dx[0], normj * dx[1], normj * dx[2]};
+
   /* signs are the same as before, since we swap i and j twice */
-  pj->primitives.gradients.rho[0] -=
-      wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-  pj->primitives.gradients.rho[1] -=
-      wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-  pj->primitives.gradients.rho[2] -=
-      wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-
-  pj->primitives.gradients.v[0][0] -=
-      wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-  pj->primitives.gradients.v[0][1] -=
-      wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-  pj->primitives.gradients.v[0][2] -=
-      wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-
-  pj->primitives.gradients.v[1][0] -=
-      wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-  pj->primitives.gradients.v[1][1] -=
-      wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-  pj->primitives.gradients.v[1][2] -=
-      wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-  pj->primitives.gradients.v[2][0] -=
-      wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-  pj->primitives.gradients.v[2][1] -=
-      wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-  pj->primitives.gradients.v[2][2] -=
-      wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-
-  pj->primitives.gradients.P[0] -=
-      wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
-  pj->primitives.gradients.P[1] -=
-      wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
-  pj->primitives.gradients.P[2] -=
-      wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+  pj->gradients.rho[0] -= dW[0] * njdx[0];
+  pj->gradients.rho[1] -= dW[0] * njdx[1];
+  pj->gradients.rho[2] -= dW[0] * njdx[2];
+
+  pj->gradients.v[0][0] -= dW[1] * njdx[0];
+  pj->gradients.v[0][1] -= dW[1] * njdx[1];
+  pj->gradients.v[0][2] -= dW[1] * njdx[2];
+  pj->gradients.v[1][0] -= dW[2] * njdx[0];
+  pj->gradients.v[1][1] -= dW[2] * njdx[1];
+  pj->gradients.v[1][2] -= dW[2] * njdx[2];
+  pj->gradients.v[2][0] -= dW[3] * njdx[0];
+  pj->gradients.v[2][1] -= dW[3] * njdx[1];
+  pj->gradients.v[2][2] -= dW[3] * njdx[2];
+
+  pj->gradients.P[0] -= dW[4] * njdx[0];
+  pj->gradients.P[1] -= dW[4] * njdx[1];
+  pj->gradients.P[2] -= dW[4] * njdx[2];
 
   hydro_slope_limit_cell_collect(pj, pi, r);
 }
@@ -169,7 +144,7 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
                                struct part *restrict pi,
                                struct part *restrict pj) {
 
-  const float r_inv = 1.f / sqrtf(r2);
+  const float r_inv = 1.0f / sqrtf(r2);
   const float r = r2 * r_inv;
 
   float wi, wi_dx;
@@ -177,41 +152,29 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
   const float xi = r * hi_inv;
   kernel_deval(xi, &wi, &wi_dx);
 
-  /* very basic gradient estimate */
-  pi->primitives.gradients.rho[0] -=
-      wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-  pi->primitives.gradients.rho[1] -=
-      wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-  pi->primitives.gradients.rho[2] -=
-      wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-
-  pi->primitives.gradients.v[0][0] -=
-      wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-  pi->primitives.gradients.v[0][1] -=
-      wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-  pi->primitives.gradients.v[0][2] -=
-      wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-
-  pi->primitives.gradients.v[1][0] -=
-      wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-  pi->primitives.gradients.v[1][1] -=
-      wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-  pi->primitives.gradients.v[1][2] -=
-      wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-
-  pi->primitives.gradients.v[2][0] -=
-      wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-  pi->primitives.gradients.v[2][1] -=
-      wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-  pi->primitives.gradients.v[2][2] -=
-      wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-
-  pi->primitives.gradients.P[0] -=
-      wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
-  pi->primitives.gradients.P[1] -=
-      wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
-  pi->primitives.gradients.P[2] -=
-      wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+  const float dW[5] = {pi->rho - pj->rho, pi->v[0] - pj->v[0],
+                       pi->v[1] - pj->v[1], pi->v[2] - pj->v[2], pi->P - pj->P};
+
+  const float normi = wi_dx * r_inv;
+  const float nidx[3] = {normi * dx[0], normi * dx[1], normi * dx[2]};
+
+  pi->gradients.rho[0] -= dW[0] * nidx[0];
+  pi->gradients.rho[1] -= dW[0] * nidx[1];
+  pi->gradients.rho[2] -= dW[0] * nidx[2];
+
+  pi->gradients.v[0][0] -= dW[1] * nidx[0];
+  pi->gradients.v[0][1] -= dW[1] * nidx[1];
+  pi->gradients.v[0][2] -= dW[1] * nidx[2];
+  pi->gradients.v[1][0] -= dW[2] * nidx[0];
+  pi->gradients.v[1][1] -= dW[2] * nidx[1];
+  pi->gradients.v[1][2] -= dW[2] * nidx[2];
+  pi->gradients.v[2][0] -= dW[3] * nidx[0];
+  pi->gradients.v[2][1] -= dW[3] * nidx[1];
+  pi->gradients.v[2][2] -= dW[3] * nidx[2];
+
+  pi->gradients.P[0] -= dW[4] * nidx[0];
+  pi->gradients.P[1] -= dW[4] * nidx[1];
+  pi->gradients.P[2] -= dW[4] * nidx[2];
 
   hydro_slope_limit_cell_collect(pi, pj, r);
 }
@@ -229,26 +192,28 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
   const float ihdimp1 = pow_dimension_plus_one(ih);
   const float volume = p->geometry.volume;
 
+  const float norm = ihdimp1 * volume;
+
   /* finalize gradients by multiplying with volume */
-  p->primitives.gradients.rho[0] *= ihdimp1 * volume;
-  p->primitives.gradients.rho[1] *= ihdimp1 * volume;
-  p->primitives.gradients.rho[2] *= ihdimp1 * volume;
+  p->gradients.rho[0] *= norm;
+  p->gradients.rho[1] *= norm;
+  p->gradients.rho[2] *= norm;
 
-  p->primitives.gradients.v[0][0] *= ihdimp1 * volume;
-  p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
-  p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
+  p->gradients.v[0][0] *= norm;
+  p->gradients.v[0][1] *= norm;
+  p->gradients.v[0][2] *= norm;
 
-  p->primitives.gradients.v[1][0] *= ihdimp1 * volume;
-  p->primitives.gradients.v[1][1] *= ihdimp1 * volume;
-  p->primitives.gradients.v[1][2] *= ihdimp1 * volume;
+  p->gradients.v[1][0] *= norm;
+  p->gradients.v[1][1] *= norm;
+  p->gradients.v[1][2] *= norm;
 
-  p->primitives.gradients.v[2][0] *= ihdimp1 * volume;
-  p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
-  p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
+  p->gradients.v[2][0] *= norm;
+  p->gradients.v[2][1] *= norm;
+  p->gradients.v[2][2] *= norm;
 
-  p->primitives.gradients.P[0] *= ihdimp1 * volume;
-  p->primitives.gradients.P[1] *= ihdimp1 * volume;
-  p->primitives.gradients.P[2] *= ihdimp1 * volume;
+  p->gradients.P[0] *= norm;
+  p->gradients.P[1] *= norm;
+  p->gradients.P[2] *= norm;
 
   hydro_slope_limit_cell(p);
 }

src/hydro/GizmoMFM/hydro_iact.h

@@ -57,7 +57,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   const float r = sqrtf(r2);
 
   /* Compute density of pi. */
-  const float hi_inv = 1.f / hi;
+  const float hi_inv = 1.0f / hi;
   const float xi = r * hi_inv;
   kernel_deval(xi, &wi, &wi_dx);
 
@@ -75,7 +75,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   pi->geometry.centroid[2] -= dx[2] * wi;
 
   /* Compute density of pj. */
-  const float hj_inv = 1.f / hj;
+  const float hj_inv = 1.0f / hj;
   const float xj = r * hj_inv;
   kernel_deval(xj, &wj, &wj_dx);
 
@@ -123,7 +123,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   /* Get r and h inverse. */
   const float r = sqrtf(r2);
 
-  const float hi_inv = 1.f / hi;
+  const float hi_inv = 1.0f / hi;
   const float xi = r * hi_inv;
   kernel_deval(xi, &wi, &wi_dx);
 
@@ -220,7 +220,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
     float r2, const float *dx, float hi, float hj, struct part *restrict pi,
     struct part *restrict pj, int mode, float a, float H) {
 
-  const float r_inv = 1.f / sqrtf(r2);
+  const float r_inv = 1.0f / sqrtf(r2);
   const float r = r2 * r_inv;
 
   /* Initialize local variables */
@@ -238,16 +238,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   const float Vi = pi->geometry.volume;
   const float Vj = pj->geometry.volume;
   float Wi[5], Wj[5];
-  Wi[0] = pi->primitives.rho;
-  Wi[1] = pi->primitives.v[0];
-  Wi[2] = pi->primitives.v[1];
-  Wi[3] = pi->primitives.v[2];
-  Wi[4] = pi->primitives.P;
-  Wj[0] = pj->primitives.rho;
-  Wj[1] = pj->primitives.v[0];
-  Wj[2] = pj->primitives.v[1];
-  Wj[3] = pj->primitives.v[2];
-  Wj[4] = pj->primitives.P;
+  Wi[0] = pi->rho;
+  Wi[1] = pi->v[0];
+  Wi[2] = pi->v[1];
+  Wi[3] = pi->v[2];
+  Wi[4] = pi->P;
+  Wj[0] = pj->rho;
+  Wj[1] = pj->v[0];
+  Wj[2] = pj->v[1];
+  Wj[3] = pj->v[2];
+  Wj[4] = pj->P;
 
   /* calculate the maximal signal velocity */
   float vmax;
@@ -258,21 +258,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   } else
     vmax = 0.0f;
 
-  float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
-               (pi->v[2] - pj->v[2]) * dx[2];
-
   /* Velocity on the axis linking the particles */
-  float dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
-                  (Wi[3] - Wj[3]) * dx[2];
-  dvdotdx = min(dvdotdx, dvdr);
+  float dvdr = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
+               (Wi[3] - Wj[3]) * dx[2];
 
-  /* We only care about this velocity for particles moving towards each others
+  /* We only care about this velocity for particles moving towards each other
    */
-  dvdotdx = min(dvdotdx, 0.f);
+  const float dvdotdx = min(dvdr, 0.0f);
 
   /* Get the signal velocity */
   /* the magical factor 3 also appears in Gadget2 */
-  vmax -= 3.f * dvdotdx * r_inv;
+  vmax -= 3.0f * dvdotdx * r_inv;
 
   /* Store the signal velocity */
   pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
@@ -280,14 +276,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 
   /* Compute kernel of pi. */
   float wi, wi_dx;
-  const float hi_inv = 1.f / hi;
+  const float hi_inv = 1.0f / hi;
   const float hi_inv_dim = pow_dimension(hi_inv);
   const float xi = r * hi_inv;
   kernel_deval(xi, &wi, &wi_dx);
 
   /* Compute kernel of pj. */
   float wj, wj_dx;
-  const float hj_inv = 1.f / hj;
+  const float hj_inv = 1.0f / hj;
   const float hj_inv_dim = pow_dimension(hj_inv);
   const float xj = r * hj_inv;
   kernel_deval(xj, &wj, &wj_dx);
@@ -298,17 +294,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   const float wi_dr = hidp1 * wi_dx;
   const float wj_dr = hjdp1 * wj_dx;
   dvdr *= r_inv;
-  if (pj->primitives.rho > 0.)
-    pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
-  if (mode == 1 && pi->primitives.rho > 0.)
-    pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
+  if (pj->rho > 0.0f)
+    pi->force.h_dt -= pj->conserved.mass * dvdr / pj->rho * wi_dr;
+  if (mode == 1 && pi->rho > 0.0f)
+    pj->force.h_dt -= pi->conserved.mass * dvdr / pi->rho * wj_dr;
 
   /* Compute (square of) area */
   /* eqn. (7) */
   float Anorm2 = 0.0f;
   float A[3];
-  if (pi->density.wcorr > const_gizmo_min_wcorr &&
-      pj->density.wcorr > const_gizmo_min_wcorr) {
+  if (pi->geometry.wcorr > const_gizmo_min_wcorr &&
+      pj->geometry.wcorr > const_gizmo_min_wcorr) {
     /* in principle, we use Vi and Vj as weights for the left and right
        contributions to the generalized surface vector.
        However, if Vi and Vj are very different (because they have very
@@ -342,10 +338,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 
   /* if the interface has no area, nothing happens and we return */
   /* continuing results in dividing by zero and NaN's... */
-  if (Anorm2 == 0.f) return;
+  if (Anorm2 == 0.0f) return;
 
   /* Compute the area */
-  const float Anorm_inv = 1. / sqrtf(Anorm2);
+  const float Anorm_inv = 1.0f / sqrtf(Anorm2);
   const float Anorm = Anorm2 * Anorm_inv;
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -407,34 +403,24 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   riemann_solve_for_middle_state_flux(Wi, Wj, n_unit, vij, totflux);
 
   /* Multiply with the interface surface area */
-  totflux[0] *= Anorm;
   totflux[1] *= Anorm;
   totflux[2] *= Anorm;
   totflux[3] *= Anorm;
   totflux[4] *= Anorm;
 
-  /* Store mass flux */
-  const float mflux_i = totflux[0];
-  pi->gravity.mflux[0] += mflux_i * dx[0];
-  pi->gravity.mflux[1] += mflux_i * dx[1];
-  pi->gravity.mflux[2] += mflux_i * dx[2];
-
   /* Update conserved variables */
+  /* We shamelessly exploit the fact that the mass flux is zero and omit all
+     terms involving it */
   /* eqn. (16) */
-  pi->conserved.flux.mass -= totflux[0];
-  pi->conserved.flux.momentum[0] -= totflux[1];
-  pi->conserved.flux.momentum[1] -= totflux[2];
-  pi->conserved.flux.momentum[2] -= totflux[3];
-  pi->conserved.flux.energy -= totflux[4];
+  pi->flux.momentum[0] -= totflux[1];
+  pi->flux.momentum[1] -= totflux[2];
+  pi->flux.momentum[2] -= totflux[3];
+  pi->flux.energy -= totflux[4];
 
 #ifndef GIZMO_TOTAL_ENERGY
-  const float ekin_i = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
-                               pi->primitives.v[1] * pi->primitives.v[1] +
-                               pi->primitives.v[2] * pi->primitives.v[2]);
-  pi->conserved.flux.energy += totflux[1] * pi->primitives.v[0];
-  pi->conserved.flux.energy += totflux[2] * pi->primitives.v[1];
-  pi->conserved.flux.energy += totflux[3] * pi->primitives.v[2];
-  pi->conserved.flux.energy -= totflux[0] * ekin_i;
+  pi->flux.energy += totflux[1] * pi->v[0];
+  pi->flux.energy += totflux[2] * pi->v[1];
+  pi->flux.energy += totflux[3] * pi->v[2];
 #endif
 
   /* Note that this used to be much more complicated in early implementations of
@@ -443,26 +429,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
    * conservation anymore and just assume the current fluxes are representative
    * for the flux over the entire time step. */
   if (mode == 1) {
-    /* Store mass flux */
-    const float mflux_j = totflux[0];
-    pj->gravity.mflux[0] -= mflux_j * dx[0];
-    pj->gravity.mflux[1] -= mflux_j * dx[1];
-    pj->gravity.mflux[2] -= mflux_j * dx[2];
-
-    pj->conserved.flux.mass += totflux[0];
-    pj->conserved.flux.momentum[0] += totflux[1];
-    pj->conserved.flux.momentum[1] += totflux[2];
-    pj->conserved.flux.momentum[2] += totflux[3];
-    pj->conserved.flux.energy += totflux[4];
+    pj->flux.momentum[0] += totflux[1];
+    pj->flux.momentum[1] += totflux[2];
+    pj->flux.momentum[2] += totflux[3];
+    pj->flux.energy += totflux[4];
 
 #ifndef GIZMO_TOTAL_ENERGY
-    const float ekin_j = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
-                                 pj->primitives.v[1] * pj->primitives.v[1] +
-                                 pj->primitives.v[2] * pj->primitives.v[2]);
-    pj->conserved.flux.energy -= totflux[1] * pj->primitives.v[0];
-    pj->conserved.flux.energy -= totflux[2] * pj->primitives.v[1];
-    pj->conserved.flux.energy -= totflux[3] * pj->primitives.v[2];
-    pj->conserved.flux.energy += totflux[0] * ekin_j;
+    pj->flux.energy -= totflux[1] * pj->v[0];
+    pj->flux.energy -= totflux[2] * pj->v[1];
+    pj->flux.energy -= totflux[3] * pj->v[2];
 #endif
   }
 }

src/hydro/GizmoMFM/hydro_io.h

@@ -63,7 +63,7 @@ INLINE static void hydro_read_particles(struct part* parts,
   list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
                                 UNIT_CONV_ACCELERATION, parts, a_hydro);
   list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
-                                UNIT_CONV_DENSITY, parts, primitives.rho);
+                                UNIT_CONV_DENSITY, parts, rho);
 }
 
 /**
@@ -203,12 +203,12 @@ INLINE static void hydro_write_particles(const struct part* parts,
                                               parts, xparts, convert_u);
   list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
                                  UNIT_CONV_NO_UNITS, parts, id);
-  list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
-                                 primitives.rho);
+  list[6] =
+      io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
   list[7] = io_make_output_field_convert_part(
       "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, xparts, convert_A);
-  list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
-                                 parts, primitives.P);
+  list[8] =
+      io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, P);
   list[9] = io_make_output_field_convert_part(
       "TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, xparts, convert_Etot);
 

src/hydro/GizmoMFM/hydro_part.h

@@ -63,31 +63,23 @@ struct part {
   /* Particle smoothing length. */
   float h;
 
-  /* The primitive hydrodynamical variables. */
-  struct {
-
-    /* Density. */
-    float rho;
-
-    /* Fluid velocity. */
-    float v[3];
+  /* Density. */
+  float rho;
 
-    /* Pressure. */
-    float P;
+  /* Pressure. */
+  float P;
 
-    /* Gradients of the primitive variables. */
+  union {
+    /* Quantities used during the volume (=density) loop. */
     struct {
 
-      /* Density gradients. */
-      float rho[3];
+      /* Derivative of particle number density. */
+      float wcount_dh;
 
-      /* Fluid velocity gradients. */
-      float v[3][3];
+      /* Particle number density. */
+      float wcount;
 
-      /* Pressure gradients. */
-      float P[3];
-
-    } gradients;
+    } density;
 
     /* Quantities needed by the slope limiter. */
     struct {
@@ -106,7 +98,56 @@ struct part {
 
     } limiter;
 
-  } primitives;
+    struct {
+      /* Fluxes. */
+      struct {
+
+        /* No mass flux, since it is always zero. */
+
+        /* Momentum flux. */
+        float momentum[3];
+
+        /* Energy flux. */
+        float energy;
+
+      } flux;
+
+      /* Variables used for timestep calculation. */
+      struct {
+
+        /* Maximum signal velocity among all the neighbours of the particle. The
+         * signal velocity encodes information about the relative fluid
+         * velocities
+         * AND particle velocities of the neighbour and this particle, as well
+         * as
+         * the sound speed of both particles. */
+        float vmax;
+
+      } timestepvars;
+
+      /* Quantities used during the force loop. */
+      struct {
+
+        /* Needed to drift the primitive variables. */
+        float h_dt;
+
+      } force;
+    };
+  };
+
+  /* Gradients of the primitive variables. */
+  struct {
+
+    /* Density gradients. */
+    float rho[3];
+
+    /* Fluid velocity gradients. */
+    float v[3][3];
+
+    /* Pressure gradients. */
+    float P[3];
+
+  } gradients;
 
   /* The conserved hydrodynamical variables. */
   struct {
@@ -120,20 +161,6 @@ struct part {
     /* Fluid thermal energy (not per unit mass!). */
     float energy;
 
-    /* Fluxes. */
-    struct {
-
-      /* Mass flux. */
-      float mass;
-
-      /* Momentum flux. */
-      float momentum[3];
-
-      /* Energy flux. */
-      float energy;
-
-    } flux;
-
   } conserved;
 
   /* Geometrical quantities used for hydro. */
@@ -149,48 +176,10 @@ struct part {
     /* Centroid of the "cell". */
     float centroid[3];
 
-  } geometry;
-
-  /* Variables used for timestep calculation. */
-  struct {
-
-    /* Maximum signal velocity among all the neighbours of the particle. The
-     * signal velocity encodes information about the relative fluid velocities
-     * AND particle velocities of the neighbour and this particle, as well as
-     * the sound speed of both particles. */
-    float vmax;
-
-  } timestepvars;
-
-  /* Quantities used during the volume (=density) loop. */
-  struct {
-
-    /* Derivative of particle number density. */
-    float wcount_dh;
-
-    /* Particle number density. */
-    float wcount;
-
     /* Correction factor for wcount. */
     float wcorr;
 
-  } density;
-
-  /* Quantities used during the force loop. */
-  struct {
-
-    /* Needed to drift the primitive variables. */
-    float h_dt;
-
-  } force;
-
-  /* Specific stuff for the gravity-hydro coupling. */
-  struct {
-
-    /* Current value of the mass flux vector. */
-    float mflux[3];
-
-  } gravity;
+  } geometry;
 
   /* Chemistry information */
   struct chemistry_part_data chemistry_data;

src/hydro/GizmoMFM/hydro_slope_limiters_cell.h

@@ -29,18 +29,18 @@
 __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
     struct part* p) {
 
-  p->primitives.limiter.rho[0] = FLT_MAX;
-  p->primitives.limiter.rho[1] = -FLT_MAX;
-  p->primitives.limiter.v[0][0] = FLT_MAX;
-  p->primitives.limiter.v[0][1] = -FLT_MAX;
-  p->primitives.limiter.v[1][0] = FLT_MAX;
-  p->primitives.limiter.v[1][1] = -FLT_MAX;
-  p->primitives.limiter.v[2][0] = FLT_MAX;
-  p->primitives.limiter.v[2][1] = -FLT_MAX;
-  p->primitives.limiter.P[0] = FLT_MAX;
-  p->primitives.limiter.P[1] = -FLT_MAX;
-
-  p->primitives.limiter.maxr = -FLT_MAX;
+  p->limiter.rho[0] = FLT_MAX;
+  p->limiter.rho[1] = -FLT_MAX;
+  p->limiter.v[0][0] = FLT_MAX;
+  p->limiter.v[0][1] = -FLT_MAX;
+  p->limiter.v[1][0] = FLT_MAX;
+  p->limiter.v[1][1] = -FLT_MAX;
+  p->limiter.v[2][0] = FLT_MAX;
+  p->limiter.v[2][1] = -FLT_MAX;
+  p->limiter.P[0] = FLT_MAX;
+  p->limiter.P[1] = -FLT_MAX;
+
+  p->limiter.maxr = -FLT_MAX;
 }
 
 /**
@@ -56,30 +56,20 @@ hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
 
   /* basic slope limiter: collect the maximal and the minimal value for the
    * primitive variables among the ngbs */
-  pi->primitives.limiter.rho[0] =
-      min(pj->primitives.rho, pi->primitives.limiter.rho[0]);
-  pi->primitives.limiter.rho[1] =
-      max(pj->primitives.rho, pi->primitives.limiter.rho[1]);
-
-  pi->primitives.limiter.v[0][0] =
-      min(pj->primitives.v[0], pi->primitives.limiter.v[0][0]);
-  pi->primitives.limiter.v[0][1] =
-      max(pj->primitives.v[0], pi->primitives.limiter.v[0][1]);
-  pi->primitives.limiter.v[1][0] =
-      min(pj->primitives.v[1], pi->primitives.limiter.v[1][0]);
-  pi->primitives.limiter.v[1][1] =
-      max(pj->primitives.v[1], pi->primitives.limiter.v[1][1]);
-  pi->primitives.limiter.v[2][0] =
-      min(pj->primitives.v[2], pi->primitives.limiter.v[2][0]);
-  pi->primitives.limiter.v[2][1] =
-      max(pj->primitives.v[2], pi->primitives.limiter.v[2][1]);
-
-  pi->primitives.limiter.P[0] =
-      min(pj->primitives.P, pi->primitives.limiter.P[0]);
-  pi->primitives.limiter.P[1] =
-      max(pj->primitives.P, pi->primitives.limiter.P[1]);
-
-  pi->primitives.limiter.maxr = max(r, pi->primitives.limiter.maxr);
+  pi->limiter.rho[0] = min(pj->rho, pi->limiter.rho[0]);
+  pi->limiter.rho[1] = max(pj->rho, pi->limiter.rho[1]);
+
+  pi->limiter.v[0][0] = min(pj->v[0], pi->limiter.v[0][0]);
+  pi->limiter.v[0][1] = max(pj->v[0], pi->limiter.v[0][1]);
+  pi->limiter.v[1][0] = min(pj->v[1], pi->limiter.v[1][0]);
+  pi->limiter.v[1][1] = max(pj->v[1], pi->limiter.v[1][1]);
+  pi->limiter.v[2][0] = min(pj->v[2], pi->limiter.v[2][0]);
+  pi->limiter.v[2][1] = max(pj->v[2], pi->limiter.v[2][1]);
+
+  pi->limiter.P[0] = min(pj->P, pi->limiter.P[0]);
+  pi->limiter.P[1] = max(pj->P, pi->limiter.P[1]);
+
+  pi->limiter.maxr = max(r, pi->limiter.maxr);
 }
 
 /**
@@ -92,94 +82,94 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
 
   float gradtrue, gradrho[3], gradv[3][3], gradP[3];
 
-  gradrho[0] = p->primitives.gradients.rho[0];
-  gradrho[1] = p->primitives.gradients.rho[1];
-  gradrho[2] = p->primitives.gradients.rho[2];
+  gradrho[0] = p->gradients.rho[0];
+  gradrho[1] = p->gradients.rho[1];
+  gradrho[2] = p->gradients.rho[2];
 
-  gradv[0][0] = p->primitives.gradients.v[0][0];
-  gradv[0][1] = p->primitives.gradients.v[0][1];
-  gradv[0][2] = p->primitives.gradients.v[0][2];
+  gradv[0][0] = p->gradients.v[0][0];
+  gradv[0][1] = p->gradients.v[0][1];
+  gradv[0][2] = p->gradients.v[0][2];
 
-  gradv[1][0] = p->primitives.gradients.v[1][0];
-  gradv[1][1] = p->primitives.gradients.v[1][1];
-  gradv[1][2] = p->primitives.gradients.v[1][2];
+  gradv[1][0] = p->gradients.v[1][0];
+  gradv[1][1] = p->gradients.v[1][1];
+  gradv[1][2] = p->gradients.v[1][2];
 
-  gradv[2][0] = p->primitives.gradients.v[2][0];
-  gradv[2][1] = p->primitives.gradients.v[2][1];
-  gradv[2][2] = p->primitives.gradients.v[2][2];
+  gradv[2][0] = p->gradients.v[2][0];
+  gradv[2][1] = p->gradients.v[2][1];
+  gradv[2][2] = p->gradients.v[2][2];
 
-  gradP[0] = p->primitives.gradients.P[0];
-  gradP[1] = p->primitives.gradients.P[1];
-  gradP[2] = p->primitives.gradients.P[2];
+  gradP[0] = p->gradients.P[0];
+  gradP[1] = p->gradients.P[1];
+  gradP[2] = p->gradients.P[2];
 
   gradtrue = sqrtf(gradrho[0] * gradrho[0] + gradrho[1] * gradrho[1] +
                    gradrho[2] * gradrho[2]);
   if (gradtrue) {
-    gradtrue *= p->primitives.limiter.maxr;
-    const float gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
-    const float gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
-    const float gradtrue_inv = 1.f / gradtrue;
+    gradtrue *= p->limiter.maxr;
+    const float gradmax = p->limiter.rho[1] - p->rho;
+    const float gradmin = p->rho - p->limiter.rho[0];
+    const float gradtrue_inv = 1.0f / gradtrue;
     const float alpha =
         min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
-    p->primitives.gradients.rho[0] *= alpha;
-    p->primitives.gradients.rho[1] *= alpha;
-    p->primitives.gradients.rho[2] *= alpha;
+    p->gradients.rho[0] *= alpha;
+    p->gradients.rho[1] *= alpha;
+    p->gradients.rho[2] *= alpha;
   }
 
   gradtrue = sqrtf(gradv[0][0] * gradv[0][0] + gradv[0][1] * gradv[0][1] +
                    gradv[0][2] * gradv[0][2]);
   if (gradtrue) {
-    gradtrue *= p->primitives.limiter.maxr;
-    const float gradmax = p->primitives.limiter.v[0][1] - p->primitives.v[0];
-    const float gradmin = p->primitives.v[0] - p->primitives.limiter.v[0][0];
-    const float gradtrue_inv = 1.f / gradtrue;
+    gradtrue *= p->limiter.maxr;
+    const float gradmax = p->limiter.v[0][1] - p->v[0];
+    const float gradmin = p->v[0] - p->limiter.v[0][0];
+    const float gradtrue_inv = 1.0f / gradtrue;
     const float alpha =
         min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
-    p->primitives.gradients.v[0][0] *= alpha;
-    p->primitives.gradients.v[0][1] *= alpha;
-    p->primitives.gradients.v[0][2] *= alpha;
+    p->gradients.v[0][0] *= alpha;
+    p->gradients.v[0][1] *= alpha;
+    p->gradients.v[0][2] *= alpha;
   }
 
   gradtrue = sqrtf(gradv[1][0] * gradv[1][0] + gradv[1][1] * gradv[1][1] +
                    gradv[1][2] * gradv[1][2]);
   if (gradtrue) {
-    gradtrue *= p->primitives.limiter.maxr;
-    const float gradmax = p->primitives.limiter.v[1][1] - p->primitives.v[1];
-    const float gradmin = p->primitives.v[1] - p->primitives.limiter.v[1][0];
-    const float gradtrue_inv = 1.f / gradtrue;
+    gradtrue *= p->limiter.maxr;
+    const float gradmax = p->limiter.v[1][1] - p->v[1];
+    const float gradmin = p->v[1] - p->limiter.v[1][0];
+    const float gradtrue_inv = 1.0f / gradtrue;
     const float alpha =
         min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
-    p->primitives.gradients.v[1][0] *= alpha;
-    p->primitives.gradients.v[1][1] *= alpha;
-    p->primitives.gradients.v[1][2] *= alpha;
+    p->gradients.v[1][0] *= alpha;
+    p->gradients.v[1][1] *= alpha;
+    p->gradients.v[1][2] *= alpha;
   }
 
   gradtrue = sqrtf(gradv[2][0] * gradv[2][0] + gradv[2][1] * gradv[2][1] +
                    gradv[2][2] * gradv[2][2]);
   if (gradtrue) {
-    gradtrue *= p->primitives.limiter.maxr;
-    const float gradmax = p->primitives.limiter.v[2][1] - p->primitives.v[2];
-    const float gradmin = p->primitives.v[2] - p->primitives.limiter.v[2][0];
-    const float gradtrue_inv = 1.f / gradtrue;
+    gradtrue *= p->limiter.maxr;
+    const float gradmax = p->limiter.v[2][1] - p->v[2];
+    const float gradmin = p->v[2] - p->limiter.v[2][0];
+    const float gradtrue_inv = 1.0f / gradtrue;
     const float alpha =
         min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
-    p->primitives.gradients.v[2][0] *= alpha;
-    p->primitives.gradients.v[2][1] *= alpha;
-    p->primitives.gradients.v[2][2] *= alpha;
+    p->gradients.v[2][0] *= alpha;
+    p->gradients.v[2][1] *= alpha;
+    p->gradients.v[2][2] *= alpha;
   }
 
   gradtrue =
       sqrtf(gradP[0] * gradP[0] + gradP[1] * gradP[1] + gradP[2] * gradP[2]);
   if (gradtrue) {
-    gradtrue *= p->primitives.limiter.maxr;
-    const float gradmax = p->primitives.limiter.P[1] - p->primitives.P;
-    const float gradmin = p->primitives.P - p->primitives.limiter.P[0];
-    const float gradtrue_inv = 1.f / gradtrue;
+    gradtrue *= p->limiter.maxr;
+    const float gradmax = p->limiter.P[1] - p->P;
+    const float gradmin = p->P - p->limiter.P[0];
+    const float gradtrue_inv = 1.0f / gradtrue;
     const float alpha =
         min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
-    p->primitives.gradients.P[0] *= alpha;
-    p->primitives.gradients.P[1] *= alpha;
-    p->primitives.gradients.P[2] *= alpha;
+    p->gradients.P[0] *= alpha;
+    p->gradients.P[1] *= alpha;
+    p->gradients.P[2] *= alpha;
   }
 }
 

src/hydro/GizmoMFM/hydro_slope_limiters_face.h

@@ -56,15 +56,17 @@ hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
 
   float phiplus, phiminus, phi_mid;
 
-  if (same_signf(phimax + delta1, phimax))
+  if (same_signf(phimax + delta1, phimax)) {
     phiplus = phimax + delta1;
-  else
+  } else {
     phiplus = phimax / (1.0f + delta1 / (fabsf(phimax) + FLT_MIN));
+  }
 
-  if (same_signf(phimin - delta1, phimin))
+  if (same_signf(phimin - delta1, phimin)) {
     phiminus = phimin - delta1;
-  else
+  } else {
     phiminus = phimin / (1.0f + delta1 / (fabsf(phimin) + FLT_MIN));
+  }
 
   if (phi_i < phi_j) {
     const float temp = min(phibar + delta2, phi_mid0);

src/hydro/GizmoMFM/hydro_velocities.h

-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_HYDRO_VELOCITIES_H
-#define SWIFT_HYDRO_VELOCITIES_H
-
-/**
- * @brief Initialize the GIZMO particle velocities before the start of the
- * actual run based on the initial value of the primitive velocity.
- *
- * @param p The particle to act upon.
- * @param xp The extended particle data to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_velocities_init(
-    struct part* restrict p, struct xpart* restrict xp) {
-
-#ifdef GIZMO_FIX_PARTICLES
-  p->v[0] = 0.f;
-  p->v[1] = 0.f;
-  p->v[2] = 0.f;
-#else
-  p->v[0] = p->primitives.v[0];
-  p->v[1] = p->primitives.v[1];
-  p->v[2] = p->primitives.v[2];
-#endif
-
-  xp->v_full[0] = p->v[0];
-  xp->v_full[1] = p->v[1];
-  xp->v_full[2] = p->v[2];
-}
-
-/**
- * @brief Set the particle velocity field that will be used to deboost fluid
- * velocities during the force loop.
- *
- * @param p The particle to act upon.
- * @param xp The extended particle data to act upon.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_velocities_prepare_force(struct part* restrict p,
-                               const struct xpart* restrict xp) {}
-
-/**
- * @brief Set the variables that will be used to update the smoothing length
- * during the drift (these will depend on the movement of the particles).
- *
- * @param p The particle to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_velocities_end_force(
-    struct part* restrict p) {
-
-#ifdef GIZMO_FIX_PARTICLES
-  /* disable the smoothing length update, since the smoothing lengths should
-     stay the same for all steps (particles don't move) */
-  p->force.h_dt = 0.0f;
-#else
-  /* Add normalization to h_dt. */
-  p->force.h_dt *= p->h * hydro_dimension_inv;
-#endif
-}
-
-/**
- * @brief Set the velocity of a GIZMO particle, based on the values of its
- * primitive variables and the geometry of its mesh-free "cell".
- *
- * @param p The particle to act upon.
- * @param xp The extended particle data to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_velocities_set(
-    struct part* restrict p, struct xpart* restrict xp) {
-
-/* We first set the particle velocity. */
-#ifdef GIZMO_FIX_PARTICLES
-
-  p->v[0] = 0.f;
-  p->v[1] = 0.f;
-  p->v[2] = 0.f;
-
-#else  // GIZMO_FIX_PARTICLES
-
-  if (p->conserved.mass > 0.f && p->primitives.rho > 0.f) {
-
-    const float inverse_mass = 1.f / p->conserved.mass;
-
-    /* Normal case: set particle velocity to fluid velocity. */
-    p->v[0] = p->conserved.momentum[0] * inverse_mass;
-    p->v[1] = p->conserved.momentum[1] * inverse_mass;
-    p->v[2] = p->conserved.momentum[2] * inverse_mass;
-
-#ifdef GIZMO_STEER_MOTION
-
-    /* Add a correction to the velocity to keep particle positions close enough
-       to
-       the centroid of their mesh-free "cell". */
-    /* The correction term below is the same one described in Springel (2010).
-     */
-    float ds[3];
-    ds[0] = p->geometry.centroid[0];
-    ds[1] = p->geometry.centroid[1];
-    ds[2] = p->geometry.centroid[2];
-    const float d = sqrtf(ds[0] * ds[0] + ds[1] * ds[1] + ds[2] * ds[2]);
-    const float R = get_radius_dimension_sphere(p->geometry.volume);
-    const float eta = 0.25f;
-    const float etaR = eta * R;
-    const float xi = 1.f;
-    const float soundspeed =
-        sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
-    /* We only apply the correction if the offset between centroid and position
-       is too large. */
-    if (d > 0.9f * etaR) {
-      float fac = xi * soundspeed / d;
-      if (d < 1.1f * etaR) {
-        fac *= 5.f * (d - 0.9f * etaR) / etaR;
-      }
-      p->v[0] -= ds[0] * fac;
-      p->v[1] -= ds[1] * fac;
-      p->v[2] -= ds[2] * fac;
-    }
-
-#endif  // GIZMO_STEER_MOTION
-  } else {
-    /* Vacuum particles have no fluid velocity. */
-    p->v[0] = 0.f;
-    p->v[1] = 0.f;
-    p->v[2] = 0.f;
-  }
-
-#endif  // GIZMO_FIX_PARTICLES
-
-  /* Now make sure all velocity variables are up to date. */
-  xp->v_full[0] = p->v[0];
-  xp->v_full[1] = p->v[1];
-  xp->v_full[2] = p->v[2];
-
-  if (p->gpart) {
-    p->gpart->v_full[0] = p->v[0];
-    p->gpart->v_full[1] = p->v[1];
-    p->gpart->v_full[2] = p->v[2];
-  }
-}
-
-#endif /* SWIFT_HYDRO_VELOCITIES_H */

src/riemann/riemann_hllc.h

@@ -33,11 +33,6 @@
 #include "riemann_checks.h"
 #include "riemann_vacuum.h"
 
-#ifndef EOS_IDEAL_GAS
-#error \
-    "The HLLC Riemann solver currently only supports and ideal gas equation of state. Either select this equation of state, or try using another Riemann solver!"
-#endif
-
 __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
     const float *WL, const float *WR, const float *n, const float *vij,
     float *totflux) {
@@ -48,19 +43,21 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
 
   /* Handle pure vacuum */
   if (!WL[0] && !WR[0]) {
-    totflux[0] = 0.f;
-    totflux[1] = 0.f;
-    totflux[2] = 0.f;
-    totflux[3] = 0.f;
-    totflux[4] = 0.f;
+    totflux[0] = 0.0f;
+    totflux[1] = 0.0f;
+    totflux[2] = 0.0f;
+    totflux[3] = 0.0f;
+    totflux[4] = 0.0f;
     return;
   }
 
   /* STEP 0: obtain velocity in interface frame */
   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]);
+  const float rhoLinv = 1.0f / WL[0];
+  const float rhoRinv = 1.0f / WR[0];
+  const float aL = sqrtf(hydro_gamma * WL[4] * rhoLinv);
+  const float aR = sqrtf(hydro_gamma * WR[4] * rhoRinv);
 
   /* Handle vacuum: vacuum does not require iteration and is always exact */
   if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) {
@@ -69,98 +66,90 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
   }
 
   /* STEP 1: pressure estimate */
-  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);
+  const float rhobar = WL[0] + WR[0];
+  const float abar = aL + aR;
+  const float pPVRS =
+      0.5f * ((WL[4] + WR[4]) - 0.25f * (uR - uL) * rhobar * abar);
+  const float pstar = max(0.0f, pPVRS);
 
   /* STEP 2: wave speed estimates
      all these speeds are along the interface normal, since uL and uR are */
-  float qL = 1.f;
-  if (pstar > WL[4] && WL[4] > 0.f) {
-    qL = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
-                         (pstar / WL[4] - 1.f));
+  float qL = 1.0f;
+  if (pstar > WL[4] && WL[4] > 0.0f) {
+    qL = sqrtf(1.0f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+                          (pstar / WL[4] - 1.0f));
   }
-  float qR = 1.f;
-  if (pstar > WR[4] && WR[4] > 0.f) {
-    qR = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
-                         (pstar / WR[4] - 1.f));
+  float qR = 1.0f;
+  if (pstar > WR[4] && WR[4] > 0.0f) {
+    qR = sqrtf(1.0f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+                          (pstar / WR[4] - 1.0f));
   }
-  const float SL = uL - aL * qL;
-  const float SR = uR + aR * qR;
+  const float SLmuL = -aL * qL;
+  const float SRmuR = 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));
+      (WR[4] - WL[4] + WL[0] * uL * SLmuL - WR[0] * uR * SRmuR) /
+      (WL[0] * SLmuL - WR[0] * SRmuR);
 
   /* STEP 3: HLLC flux in a frame moving with the interface velocity */
-  if (Sstar >= 0.f) {
+  if (Sstar >= 0.0f) {
+    const float rhoLuL = WL[0] * uL;
+    const float v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3];
+    const float eL =
+        WL[4] * rhoLinv * hydro_one_over_gamma_minus_one + 0.5f * v2;
+    const float SL = SLmuL + uL;
+
     /* flux FL */
-    totflux[0] = WL[0] * uL;
+    totflux[0] = rhoLuL;
     /* these are the actual correct fluxes in the boosted lab frame
        (not rotated to interface frame) */
-    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];
-    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.f) {
-
-      float UstarL[5];
-
-      /* add flux FstarL */
-      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 */
-      UstarL[1] = WL[1] + (Sstar - uL) * n[0];
-      UstarL[2] = WL[2] + (Sstar - uL) * n[1];
-      UstarL[3] = WL[3] + (Sstar - uL) * n[2];
-      UstarL[4] = eL + (Sstar - uL) * (Sstar + WL[4] / (WL[0] * (SL - uL)));
-      UstarL[0] *= WL[0] * (SL - uL) / (SL - Sstar);
-      UstarL[1] *= WL[0] * (SL - uL) / (SL - Sstar);
-      UstarL[2] *= WL[0] * (SL - uL) / (SL - Sstar);
-      UstarL[3] *= WL[0] * (SL - uL) / (SL - Sstar);
-      UstarL[4] *= WL[0] * (SL - uL) / (SL - Sstar);
-      totflux[0] += SL * (UstarL[0] - WL[0]);
-      totflux[1] += SL * (UstarL[1] - WL[0] * WL[1]);
-      totflux[2] += SL * (UstarL[2] - WL[0] * WL[2]);
-      totflux[3] += SL * (UstarL[3] - WL[0] * WL[3]);
-      totflux[4] += SL * (UstarL[4] - WL[0] * eL);
+    totflux[1] = rhoLuL * WL[1] + WL[4] * n[0];
+    totflux[2] = rhoLuL * WL[2] + WL[4] * n[1];
+    totflux[3] = rhoLuL * WL[3] + WL[4] * n[2];
+    totflux[4] = rhoLuL * eL + WL[4] * uL;
+
+    if (SL < 0.0f) {
+
+      const float starfac = SLmuL / (SL - Sstar) - 1.0f;
+      const float rhoLSL = WL[0] * SL;
+      const float SstarmuL = Sstar - uL;
+      const float rhoLSLstarfac = rhoLSL * starfac;
+      const float rhoLSLSstarmuL = rhoLSL * SstarmuL;
+
+      totflux[0] += rhoLSLstarfac;
+      totflux[1] += rhoLSLstarfac * WL[1] + rhoLSLSstarmuL * n[0];
+      totflux[2] += rhoLSLstarfac * WL[2] + rhoLSLSstarmuL * n[1];
+      totflux[3] += rhoLSLstarfac * WL[3] + rhoLSLSstarmuL * n[2];
+      totflux[4] += rhoLSLstarfac * eL +
+                    rhoLSLSstarmuL * (Sstar + WL[4] / (WL[0] * SLmuL));
     }
   } else {
-    /* flux FR */
-    totflux[0] = WR[0] * uR;
-    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];
+    const float rhoRuR = WR[0] * uR;
     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.f) {
+    const float eR =
+        WR[4] * rhoRinv * hydro_one_over_gamma_minus_one + 0.5f * v2;
+    const float SR = SRmuR + uR;
 
-      float UstarR[5];
-
-      /* add flux FstarR */
-      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];
-      UstarR[4] = eR + (Sstar - uR) * (Sstar + WR[4] / (WR[0] * (SR - uR)));
-      UstarR[0] *= WR[0] * (SR - uR) / (SR - Sstar);
-      UstarR[1] *= WR[0] * (SR - uR) / (SR - Sstar);
-      UstarR[2] *= WR[0] * (SR - uR) / (SR - Sstar);
-      UstarR[3] *= WR[0] * (SR - uR) / (SR - Sstar);
-      UstarR[4] *= WR[0] * (SR - uR) / (SR - Sstar);
-      totflux[0] += SR * (UstarR[0] - WR[0]);
-      totflux[1] += SR * (UstarR[1] - WR[0] * WR[1]);
-      totflux[2] += SR * (UstarR[2] - WR[0] * WR[2]);
-      totflux[3] += SR * (UstarR[3] - WR[0] * WR[3]);
-      totflux[4] += SR * (UstarR[4] - WR[0] * eR);
+    /* flux FR */
+    totflux[0] = rhoRuR;
+    totflux[1] = rhoRuR * WR[1] + WR[4] * n[0];
+    totflux[2] = rhoRuR * WR[2] + WR[4] * n[1];
+    totflux[3] = rhoRuR * WR[3] + WR[4] * n[2];
+    totflux[4] = rhoRuR * eR + WR[4] * uR;
+
+    if (SR > 0.0f) {
+
+      const float starfac = SRmuR / (SR - Sstar) - 1.0f;
+      const float rhoRSR = WR[0] * SR;
+      const float SstarmuR = Sstar - uR;
+      const float rhoRSRstarfac = rhoRSR * starfac;
+      const float rhoRSRSstarmuR = rhoRSR * SstarmuR;
+
+      totflux[0] += rhoRSRstarfac;
+      totflux[1] += rhoRSRstarfac * WR[1] + rhoRSRSstarmuR * n[0];
+      totflux[2] += rhoRSRstarfac * WR[2] + rhoRSRSstarmuR * n[1];
+      totflux[3] += rhoRSRstarfac * WR[3] + rhoRSRSstarmuR * n[2];
+      totflux[4] += rhoRSRstarfac * eR +
+                    rhoRSRSstarmuR * (Sstar + WR[4] / (WR[0] * SRmuR));
     }
   }
 
@@ -194,11 +183,11 @@ riemann_solve_for_middle_state_flux(const float *WL, const float *WR,
 
   /* Handle pure vacuum */
   if (!WL[0] && !WR[0]) {
-    totflux[0] = 0.f;
-    totflux[1] = 0.f;
-    totflux[2] = 0.f;
-    totflux[3] = 0.f;
-    totflux[4] = 0.f;
+    totflux[0] = 0.0f;
+    totflux[1] = 0.0f;
+    totflux[2] = 0.0f;
+    totflux[3] = 0.0f;
+    totflux[4] = 0.0f;
     return;
   }
 
@@ -210,37 +199,38 @@ riemann_solve_for_middle_state_flux(const float *WL, const float *WR,
 
   /* Handle vacuum: vacuum does not require iteration and is always exact */
   if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) {
-    totflux[0] = 0.f;
-    totflux[1] = 0.f;
-    totflux[2] = 0.f;
-    totflux[3] = 0.f;
-    totflux[4] = 0.f;
+    totflux[0] = 0.0f;
+    totflux[1] = 0.0f;
+    totflux[2] = 0.0f;
+    totflux[3] = 0.0f;
+    totflux[4] = 0.0f;
     return;
   }
 
   /* STEP 1: pressure estimate */
-  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 rhobar = WL[0] + WR[0];
+  const float abar = aL + aR;
+  const float pPVRS =
+      0.5f * ((WL[4] + WR[4]) - 0.25f * (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 */
-  float qL = 1.f;
-  if (pstar > WL[4] && WL[4] > 0.f) {
-    qL = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
-                         (pstar / WL[4] - 1.f));
+  float qL = 1.0f;
+  if (pstar > WL[4] && WL[4] > 0.0f) {
+    qL = sqrtf(1.0f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+                          (pstar / WL[4] - 1.0f));
   }
-  float qR = 1.f;
-  if (pstar > WR[4] && WR[4] > 0.f) {
-    qR = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
-                         (pstar / WR[4] - 1.f));
+  float qR = 1.0f;
+  if (pstar > WR[4] && WR[4] > 0.0f) {
+    qR = sqrtf(1.0f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+                          (pstar / WR[4] - 1.0f));
   }
-  const float SL = uL - aL * qL;
-  const float SR = uR + aR * qR;
+  const float SLmuL = -aL * qL;
+  const float SRmuR = 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));
+      (WR[4] - WL[4] + WL[0] * uL * SLmuL - WR[0] * uR * SRmuR) /
+      (WL[0] * SLmuL - WR[0] * SRmuR);
 
   totflux[0] = 0.0f;
   totflux[1] = pstar * n[0];