Implemented the changes in ls that prevent gizmo crashes in cases where the...

Implemented the changes in ls that prevent gizmo crashes in cases where the mesh is seriously distorted

src/Makefile.am

@@ -86,6 +86,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
                  hydro/Gizmo/hydro_slope_limiters_cell.h \
                  hydro/Gizmo/hydro_slope_limiters_face.h \
                  hydro/Gizmo/hydro_slope_limiters.h \
+                 hydro/Gizmo/hydro_unphysical.h \
+                 hydro/Gizmo/hydro_velocities.h \
                  hydro/Shadowswift/hydro_debug.h \
                  hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
                  hydro/Shadowswift/hydro_iact.h \

src/const.h

@@ -52,8 +52,43 @@
 /* Options to control the movement of particles for GIZMO_SPH. */
 /* This option disables particle movement */
 //#define GIZMO_FIX_PARTICLES
+/* Try to keep cells regular by adding a correction velocity. */
+#define GIZMO_STEER_MOTION
 //#define GIZMO_TOTAL_ENERGY
 
+/* Options to control handling of unphysical values (GIZMO_SPH only). */
+/* In GIZMO, mass and energy (and hence density and pressure) can in principle
+   become negative, which will cause unwanted behaviour that can make the code
+   crash.
+   If no options are selected below, we assume (and pray) that this will not
+   happen, and add no restrictions to how these variables are treated. */
+/* Check for unphysical values and crash if they occur. */
+//#define GIZMO_UNPHYSICAL_ERROR
+/* Check for unphysical values and reset them to safe values. */
+#define GIZMO_UNPHYSICAL_RESCUE
+/* Show a warning message if an unphysical value was reset (only works if
+   GIZMO_UNPHYSICAL_RESCUE is also selected). */
+#define GIZMO_UNPHYSICAL_WARNING
+
+/* Parameters that control how GIZMO handles pathological particle
+   configurations. */
+/* Show a warning message if a pathological configuration has been detected. */
+#define GIZMO_PATHOLOGICAL_WARNING
+/* Crash if a pathological configuration has been detected. */
+//#define GIZMO_PATHOLOGICAL_ERROR
+/* Maximum allowed gradient matrix condition number. If the condition number of
+   the gradient matrix (defined in equation C1 in Hopkins, 2015) is larger than
+   this value, we artificially increase the number of neighbours to get a more
+   homogeneous sampling. */
+#define const_gizmo_max_condition_number 100.0f
+/* Correction factor applied to the particle wcount to force more neighbours if
+   the condition number is too large. */
+#define const_gizmo_w_correction_factor 0.9f
+/* Lower limit on the wcount correction factor. If the condition number is still
+   too high after this wcount correction has been applied, we give up on the
+   gradient matrix and use SPH gradients instead. */
+#define const_gizmo_min_wcorr 0.5f
+
 /* Types of gradients to use for SHADOWFAX_SPH */
 /* If no option is chosen, no gradients are used (first order scheme) */
 #define SHADOWFAX_GRADIENTS

src/hydro/Gizmo/hydro.h

@@ -2,6 +2,7 @@
 /*******************************************************************************
  * This file is part of SWIFT.
  * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *               2016, 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
@@ -24,9 +25,13 @@
 #include "equation_of_state.h"
 #include "hydro_gradients.h"
 #include "hydro_space.h"
+#include "hydro_unphysical.h"
+#include "hydro_velocities.h"
 #include "minmax.h"
 #include "riemann.h"
 
+//#define GIZMO_LLOYD_ITERATION
+
 /**
  * @brief Computes the hydro time-step of a given particle
  *
@@ -40,6 +45,10 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
 
   const float CFL_condition = hydro_properties->CFL_condition;
 
+#ifdef GIZMO_LLOYD_ITERATION
+  return CFL_condition;
+#endif
+
   if (p->timestepvars.vmax == 0.) {
     /* vmax can be zero in vacuum cells that only have vacuum neighbours */
     /* in this case, the time step should be limited by the maximally
@@ -47,7 +56,9 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
        the time step to a very large value */
     return FLT_MAX;
   } else {
-    return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
+    const float psize = powf(p->geometry.volume / hydro_dimension_unit_sphere,
+                             hydro_dimension_inv);
+    return 2. * CFL_condition * psize / fabsf(p->timestepvars.vmax);
   }
 }
 
@@ -128,16 +139,27 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
                                  p->conserved.momentum[2] * p->primitives.v[2]);
 #endif
 
-#if defined(GIZMO_FIX_PARTICLES)
-  /* make sure the particles are initially at rest */
+#ifdef GIZMO_LLOYD_ITERATION
+  /* overwrite all variables to make sure they have safe values */
+  p->primitives.rho = 1.;
+  p->primitives.v[0] = 0.;
+  p->primitives.v[1] = 0.;
+  p->primitives.v[2] = 0.;
+  p->primitives.P = 1.;
+
+  p->conserved.mass = 1.;
+  p->conserved.momentum[0] = 0.;
+  p->conserved.momentum[1] = 0.;
+  p->conserved.momentum[2] = 0.;
+  p->conserved.energy = 1.;
+
   p->v[0] = 0.;
   p->v[1] = 0.;
   p->v[2] = 0.;
 #endif
 
-  xp->v_full[0] = p->v[0];
-  xp->v_full[1] = p->v[1];
-  xp->v_full[2] = p->v[2];
+  /* initialize the particle velocity based on the primitive fluid velocity */
+  hydro_velocities_init(p, xp);
 
   /* we cannot initialize wcorr in init_part, as init_part gets called every
      time the density loop is repeated, and the whole point of storing wcorr
@@ -169,6 +191,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   p->geometry.matrix_E[2][0] = 0.0f;
   p->geometry.matrix_E[2][1] = 0.0f;
   p->geometry.matrix_E[2][2] = 0.0f;
+  p->geometry.centroid[0] = 0.0f;
+  p->geometry.centroid[1] = 0.0f;
+  p->geometry.centroid[2] = 0.0f;
   p->geometry.Atot = 0.0f;
 
   /* Set the active flag to active. */
@@ -226,6 +251,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   p->geometry.matrix_E[2][1] = ihdim * p->geometry.matrix_E[2][1];
   p->geometry.matrix_E[2][2] = ihdim * p->geometry.matrix_E[2][2];
 
+  p->geometry.centroid[0] *= kernel_norm;
+  p->geometry.centroid[1] *= kernel_norm;
+  p->geometry.centroid[2] *= kernel_norm;
+
+  p->geometry.centroid[0] /= p->density.wcount;
+  p->geometry.centroid[1] /= p->density.wcount;
+  p->geometry.centroid[2] /= p->density.wcount;
+
   /* Check the condition number to see if we have a stable geometry. */
   float condition_number_E = 0.0f;
   int i, j;
@@ -249,12 +282,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   float condition_number =
       hydro_dimension_inv * sqrtf(condition_number_E * condition_number_Einv);
 
-  if (condition_number > 100.0f) {
-    //    error("Condition number larger than 100!");
-    //    message("Condition number too large: %g (p->id: %llu)!",
-    //    condition_number, p->id);
+  if (condition_number > const_gizmo_max_condition_number &&
+      p->density.wcorr > const_gizmo_min_wcorr) {
+#ifdef GIZMO_PATHOLOGICAL_ERROR
+    error("Condition number larger than %g (%g)!",
+          const_gizmo_max_condition_number, condition_number);
+#endif
+#ifdef GIZMO_PATHOLOGICAL_WARNING
+    message("Condition number too large: %g (> %g, p->id: %llu)!",
+            condition_number, const_gizmo_max_condition_number, p->id);
+#endif
     /* add a correction to the number of neighbours for this particle */
-    p->density.wcorr *= 0.75;
+    p->density.wcorr *= const_gizmo_w_correction_factor;
   }
 
   hydro_gradients_init(p);
@@ -264,8 +303,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   const float m = p->conserved.mass;
 
 #ifdef SWIFT_DEBUG_CHECKS
-  if (m == 0.) {
-    error("Mass is 0!");
+  if (m < 0.) {
+    error("Mass is negative!");
   }
 
   if (volume == 0.) {
@@ -278,15 +317,20 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   momentum[1] = p->conserved.momentum[1];
   momentum[2] = p->conserved.momentum[2];
   p->primitives.rho = m / volume;
-  p->primitives.v[0] = momentum[0] / m;
-  p->primitives.v[1] = momentum[1] / m;
-  p->primitives.v[2] = momentum[2] / m;
+  if (m == 0.) {
+    p->primitives.v[0] = 0.;
+    p->primitives.v[1] = 0.;
+    p->primitives.v[2] = 0.;
+  } else {
+    p->primitives.v[0] = momentum[0] / m;
+    p->primitives.v[1] = momentum[1] / m;
+    p->primitives.v[2] = momentum[2] / m;
+  }
 
 #ifdef EOS_ISOTHERMAL_GAS
   /* although the pressure is not formally used anywhere if an isothermal eos
      has been selected, we still make sure it is set to the correct value */
-  p->primitives.P = const_isothermal_soundspeed * const_isothermal_soundspeed *
-                    p->primitives.rho;
+  p->primitives.P = gas_pressure_from_internal_energy(p->primitives.rho, 0.);
 #else
 
   float energy = p->conserved.energy;
@@ -304,12 +348,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 #endif
 
   /* sanity checks */
-  /* it would probably be safer to throw a warning if netive densities or
-     pressures occur */
-  if (p->primitives.rho < 0.0f || p->primitives.P < 0.0f) {
-    p->primitives.rho = 0.0f;
-    p->primitives.P = 0.0f;
-  }
+  gizmo_check_physical_quantity("density", p->primitives.rho);
+  gizmo_check_physical_quantity("pressure", p->primitives.P);
+
+#ifdef GIZMO_LLOYD_ITERATION
+  /* overwrite primitive variables to make sure they still have safe values */
+  p->primitives.rho = 1.;
+  p->primitives.v[0] = 0.;
+  p->primitives.v[1] = 0.;
+  p->primitives.v[2] = 0.;
+  p->primitives.P = 1.;
+#endif
 
   /* Add a correction factor to wcount (to force a neighbour number increase if
      the geometry matrix is close to singular) */
@@ -330,8 +379,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
  *
  * @param p The particle to act upon.
  * @param xp The extended particle data to act upon.
- * @param ti_current Current integer time.
- * @param timeBase Conversion factor between integer time and physical time.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part* restrict p, struct xpart* restrict xp) {
@@ -340,10 +387,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   p->timestepvars.vmax = 0.0f;
 
   /* Set the actual velocity of the particle */
-  /* if GIZMO_FIX_PARTICLES has been selected, v_full will always be zero */
-  p->force.v_full[0] = xp->v_full[0];
-  p->force.v_full[1] = xp->v_full[1];
-  p->force.v_full[2] = xp->v_full[2];
+  hydro_velocities_prepare_force(p, xp);
 }
 
 /**
@@ -364,6 +408,11 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
   p->gravity.mflux[0] = 0.0f;
   p->gravity.mflux[1] = 0.0f;
   p->gravity.mflux[2] = 0.0f;
+
+#ifdef GIZMO_LLOYD_ITERATION
+  /* reset the gradients to zero, as we don't want them */
+  hydro_gradients_init(p);
+#endif
 }
 
 /**
@@ -422,6 +471,10 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     struct part* p, struct xpart* xp, float dt) {
 
+#ifdef GIZMO_LLOYD_ITERATION
+  return;
+#endif
+
   const float h_inv = 1.0f / p->h;
 
   /* Predict smoothing length */
@@ -432,8 +485,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   else
     h_corr = expf(w1);
 
-  /* Limit the smoothing length correction. */
-  if (h_corr < 2.0f) {
+  /* Limit the smoothing length correction (and make sure it is always
+     positive). */
+  if (h_corr < 2.0f && h_corr > 0.) {
     p->h *= h_corr;
   }
 
@@ -483,22 +537,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
 
   /* set the variables that are used to drift the primitive variables */
 
-  /* Add normalization to h_dt. */
-  p->force.h_dt *= p->h * hydro_dimension_inv;
-
-  if (p->force.dt) {
+  if (p->force.dt > 0.) {
     p->du_dt = p->conserved.flux.energy / p->force.dt;
   } else {
     p->du_dt = 0.0f;
   }
 
-#if defined(GIZMO_FIX_PARTICLES)
-  p->du_dt = 0.0f;
-
-  /* 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;
-#endif
+  hydro_velocities_end_force(p);
 }
 
 /**
@@ -527,7 +572,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   p->conserved.energy += p->conserved.flux.energy;
 #endif
 
+  gizmo_check_physical_quantity("mass", p->conserved.mass);
+  gizmo_check_physical_quantity("energy", p->conserved.energy);
+
 #ifdef SWIFT_DEBUG_CHECKS
+  /* Note that this check will only have effect if no GIZMO_UNPHYSICAL option
+     was selected. */
   if (p->conserved.mass < 0.) {
     error(
         "Negative mass after conserved variables update (mass: %g, dmass: %g)!",
@@ -535,7 +585,10 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   }
 
   if (p->conserved.energy < 0.) {
-    error("Negative energy after conserved variables update!");
+    error(
+        "Negative energy after conserved variables update (energy: %g, "
+        "denergy: %g)!",
+        p->conserved.energy, p->conserved.flux.energy);
   }
 #endif
 
@@ -549,7 +602,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     a_grav[2] = p->gpart->a_grav[2];
 
     /* Store the gravitational acceleration for later use. */
-    /* This is currently only used for output purposes. */
+    /* This is used for the prediction step. */
     p->gravity.old_a[0] = a_grav[0];
     p->gravity.old_a[1] = a_grav[1];
     p->gravity.old_a[2] = a_grav[2];
@@ -564,7 +617,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
     p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
 
-#if !defined(EOS_ISOTHERMAL_GAS) && defined(GIZMO_TOTAL_ENERGY)
+#if !defined(EOS_ISOTHERMAL_GAS)
     /* This part still needs to be tested! */
     p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] +
                                  p->conserved.momentum[1] * a_grav[1] +
@@ -585,45 +638,25 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   p->conserved.flux.momentum[2] = 0.0f;
   p->conserved.flux.energy = 0.0f;
 
-#if defined(GIZMO_FIX_PARTICLES)
-  xp->v_full[0] = 0.;
-  xp->v_full[1] = 0.;
-  xp->v_full[2] = 0.;
-
-  p->v[0] = 0.;
-  p->v[1] = 0.;
-  p->v[2] = 0.;
-
-  if (p->gpart) {
-    p->gpart->v_full[0] = 0.;
-    p->gpart->v_full[1] = 0.;
-    p->gpart->v_full[2] = 0.;
-  }
-#else
-  /* Set particle movement */
-  if (p->conserved.mass > 0.) {
-    xp->v_full[0] = p->conserved.momentum[0] / p->conserved.mass;
-    xp->v_full[1] = p->conserved.momentum[1] / p->conserved.mass;
-    xp->v_full[2] = p->conserved.momentum[2] / p->conserved.mass;
-  } else {
-    /* vacuum particles don't move */
-    xp->v_full[0] = 0.;
-    xp->v_full[1] = 0.;
-    xp->v_full[2] = 0.;
-  }
+  hydro_velocities_set(p, xp);
+
+#ifdef GIZMO_LLOYD_ITERATION
+  /* reset conserved variables to safe values */
+  p->conserved.mass = 1.;
+  p->conserved.momentum[0] = 0.;
+  p->conserved.momentum[1] = 0.;
+  p->conserved.momentum[2] = 0.;
+  p->conserved.energy = 1.;
+
+  /* set the particle velocities to the Lloyd velocities */
+  /* note that centroid is the relative position of the centroid w.r.t. the
+     particle position (position - centroid) */
+  xp->v_full[0] = -p->geometry.centroid[0] / p->force.dt;
+  xp->v_full[1] = -p->geometry.centroid[1] / p->force.dt;
+  xp->v_full[2] = -p->geometry.centroid[2] / p->force.dt;
   p->v[0] = xp->v_full[0];
   p->v[1] = xp->v_full[1];
   p->v[2] = xp->v_full[2];
-
-  /* Update gpart! */
-  /* This is essential, as the gpart drift is done independently from the part
-     drift, and we don't want the gpart and the part to have different
-     positions! */
-  if (p->gpart) {
-    p->gpart->v_full[0] = xp->v_full[0];
-    p->gpart->v_full[1] = xp->v_full[1];
-    p->gpart->v_full[2] = xp->v_full[2];
-  }
 #endif
 
   /* reset wcorr */

src/hydro/Gizmo/hydro_gradients.h

@@ -22,6 +22,7 @@
 #define SWIFT_HYDRO_GRADIENTS_H
 
 #include "hydro_slope_limiters.h"
+#include "hydro_unphysical.h"
 #include "riemann.h"
 
 #if defined(GRADIENTS_SPH)
@@ -98,6 +99,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
   float xij_j[3];
   int k;
   float xfac;
+  float a_grav_i[3], a_grav_j[3];
 
   /* perform gradient reconstruction in space and time */
   /* space */
@@ -139,37 +141,38 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
            pj->primitives.gradients.P[1] * xij_j[1] +
            pj->primitives.gradients.P[2] * xij_j[2];
 
+  a_grav_i[0] = pi->gravity.old_a[0];
+  a_grav_i[1] = pi->gravity.old_a[1];
+  a_grav_i[2] = pi->gravity.old_a[2];
+
+  a_grav_i[0] += pi->gravity.grad_a[0][0] * xij_i[0] +
+                 pi->gravity.grad_a[0][1] * xij_i[1] +
+                 pi->gravity.grad_a[0][2] * xij_i[2];
+  a_grav_i[1] += pi->gravity.grad_a[1][0] * xij_i[0] +
+                 pi->gravity.grad_a[1][1] * xij_i[1] +
+                 pi->gravity.grad_a[1][2] * xij_i[2];
+  a_grav_i[2] += pi->gravity.grad_a[2][0] * xij_i[0] +
+                 pi->gravity.grad_a[2][1] * xij_i[1] +
+                 pi->gravity.grad_a[2][2] * xij_i[2];
+
+  a_grav_j[0] = pj->gravity.old_a[0];
+  a_grav_j[1] = pj->gravity.old_a[1];
+  a_grav_j[2] = pj->gravity.old_a[2];
+
+  a_grav_j[0] += pj->gravity.grad_a[0][0] * xij_j[0] +
+                 pj->gravity.grad_a[0][1] * xij_j[1] +
+                 pj->gravity.grad_a[0][2] * xij_j[2];
+  a_grav_j[1] += pj->gravity.grad_a[1][0] * xij_j[0] +
+                 pj->gravity.grad_a[1][1] * xij_j[1] +
+                 pj->gravity.grad_a[1][2] * xij_j[2];
+  a_grav_j[2] += pj->gravity.grad_a[2][0] * xij_j[0] +
+                 pj->gravity.grad_a[2][1] * xij_j[1] +
+                 pj->gravity.grad_a[2][2] * xij_j[2];
+
   hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
 
   /* time */
   if (Wi[0] > 0.0f) {
-#ifdef EOS_ISOTHERMAL_GAS
-    dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
-                             Wi[2] * pi->primitives.gradients.rho[1] +
-                             Wi[3] * pi->primitives.gradients.rho[2] +
-                             Wi[0] * (pi->primitives.gradients.v[0][0] +
-                                      pi->primitives.gradients.v[1][1] +
-                                      pi->primitives.gradients.v[2][2]));
-    dWi[1] -= 0.5 * mindt *
-              (Wi[1] * pi->primitives.gradients.v[0][0] +
-               Wi[2] * pi->primitives.gradients.v[0][1] +
-               Wi[3] * pi->primitives.gradients.v[0][2] +
-               const_isothermal_soundspeed * const_isothermal_soundspeed *
-                   pi->primitives.gradients.rho[0] / Wi[0]);
-    dWi[2] -= 0.5 * mindt *
-              (Wi[1] * pi->primitives.gradients.v[1][0] +
-               Wi[2] * pi->primitives.gradients.v[1][1] +
-               Wi[3] * pi->primitives.gradients.v[1][2] +
-               const_isothermal_soundspeed * const_isothermal_soundspeed *
-                   pi->primitives.gradients.rho[1] / Wi[0]);
-    dWi[3] -= 0.5 * mindt *
-              (Wi[1] * pi->primitives.gradients.v[2][0] +
-               Wi[2] * pi->primitives.gradients.v[2][1] +
-               Wi[3] * pi->primitives.gradients.v[2][2] +
-               const_isothermal_soundspeed * const_isothermal_soundspeed *
-                   pi->primitives.gradients.rho[2] / Wi[0]);
-/* we don't care about P in this case */
-#else
     dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
                              Wi[2] * pi->primitives.gradients.rho[1] +
                              Wi[3] * pi->primitives.gradients.rho[2] +
@@ -195,36 +198,13 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
                hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
                                       pi->primitives.gradients.v[1][1] +
                                       pi->primitives.gradients.v[2][2]));
-#endif
+
+    dWi[1] += 0.5 * mindt * a_grav_i[0];
+    dWi[2] += 0.5 * mindt * a_grav_i[1];
+    dWi[3] += 0.5 * mindt * a_grav_i[2];
   }
 
   if (Wj[0] > 0.0f) {
-#ifdef EOS_ISOTHERMAL_GAS
-    dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
-                             Wj[2] * pj->primitives.gradients.rho[1] +
-                             Wj[3] * pj->primitives.gradients.rho[2] +
-                             Wj[0] * (pj->primitives.gradients.v[0][0] +
-                                      pj->primitives.gradients.v[1][1] +
-                                      pj->primitives.gradients.v[2][2]));
-    dWj[1] -= 0.5 * mindt *
-              (Wj[1] * pj->primitives.gradients.v[0][0] +
-               Wj[2] * pj->primitives.gradients.v[0][1] +
-               Wj[3] * pj->primitives.gradients.v[0][2] +
-               const_isothermal_soundspeed * const_isothermal_soundspeed *
-                   pj->primitives.gradients.rho[0] / Wj[0]);
-    dWj[2] -= 0.5 * mindt *
-              (Wj[1] * pj->primitives.gradients.v[1][0] +
-               Wj[2] * pj->primitives.gradients.v[1][1] +
-               Wj[3] * pj->primitives.gradients.v[1][2] +
-               const_isothermal_soundspeed * const_isothermal_soundspeed *
-                   pj->primitives.gradients.rho[1] / Wj[0]);
-    dWj[3] -= 0.5 * mindt *
-              (Wj[1] * pj->primitives.gradients.v[2][0] +
-               Wj[2] * pj->primitives.gradients.v[2][1] +
-               Wj[3] * pj->primitives.gradients.v[2][2] +
-               const_isothermal_soundspeed * const_isothermal_soundspeed *
-                   pj->primitives.gradients.rho[2] / Wj[0]);
-#else
     dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
                              Wj[2] * pj->primitives.gradients.rho[1] +
                              Wj[3] * pj->primitives.gradients.rho[2] +
@@ -250,36 +230,28 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
                hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
                                       pj->primitives.gradients.v[1][1] +
                                       pj->primitives.gradients.v[2][2]));
-#endif
-  }
 
-  if (-dWi[0] > Wi[0]) {
-    Wi[0] = 0.0f;
-  } else {
-    Wi[0] += dWi[0];
+    dWj[1] += 0.5 * mindt * a_grav_j[0];
+    dWj[2] += 0.5 * mindt * a_grav_j[1];
+    dWj[3] += 0.5 * mindt * a_grav_j[2];
   }
+
+  Wi[0] += dWi[0];
   Wi[1] += dWi[1];
   Wi[2] += dWi[2];
   Wi[3] += dWi[3];
-  if (-dWi[4] > Wi[4]) {
-    Wi[4] = 0.0f;
-  } else {
-    Wi[4] += dWi[4];
-  }
+  Wi[4] += dWi[4];
 
-  if (-dWj[0] > Wj[0]) {
-    Wj[0] = 0.0f;
-  } else {
-    Wj[0] += dWj[0];
-  }
+  Wj[0] += dWj[0];
   Wj[1] += dWj[1];
   Wj[2] += dWj[2];
   Wj[3] += dWj[3];
-  if (-dWj[4] > Wj[4]) {
-    Wj[4] = 0.0f;
-  } else {
-    Wj[4] += dWj[4];
-  }
+  Wj[4] += dWj[4];
+
+  gizmo_check_physical_quantity("density", Wi[0]);
+  gizmo_check_physical_quantity("pressure", Wi[4]);
+  gizmo_check_physical_quantity("density", Wj[0]);
+  gizmo_check_physical_quantity("pressure", Wj[4]);
 }
 
 #endif  // SWIFT_HYDRO_GRADIENTS_H

src/hydro/Gizmo/hydro_iact.h

@@ -23,6 +23,8 @@
 #include "hydro_gradients.h"
 #include "riemann.h"
 
+#define GIZMO_VOLUME_CORRECTION
+
 /**
  * @brief Calculate the volume interaction between particle i and particle j
  *
@@ -62,6 +64,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   for (k = 0; k < 3; k++)
     for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
 
+  pi->geometry.centroid[0] -= dx[0] * wi;
+  pi->geometry.centroid[1] -= dx[1] * wi;
+  pi->geometry.centroid[2] -= dx[2] * wi;
+
   /* Compute density of pj. */
   h_inv = 1.0 / hj;
   xj = r * h_inv;
@@ -74,6 +80,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   pj->geometry.volume += wj;
   for (k = 0; k < 3; k++)
     for (l = 0; l < 3; l++) pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
+
+  pj->geometry.centroid[0] += dx[0] * wj;
+  pj->geometry.centroid[1] += dx[1] * wj;
+  pj->geometry.centroid[2] += dx[2] * wj;
 }
 
 /**
@@ -117,6 +127,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   pi->geometry.volume += wi;
   for (k = 0; k < 3; k++)
     for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+
+  pi->geometry.centroid[0] -= dx[0] * wi;
+  pi->geometry.centroid[1] -= dx[1] * wi;
+  pi->geometry.centroid[2] -= dx[2] * wi;
 }
 
 /**
@@ -325,14 +339,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 
   /* calculate the maximal signal velocity */
   if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
-#ifdef EOS_ISOTHERMAL_GAS
-    /* we use a value that is slightly higher than necessary, since the correct
-       value does not always work */
-    vmax = 2.5 * const_isothermal_soundspeed;
-#else
     vmax =
         sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
-#endif
   } else {
     vmax = 0.0f;
   }
@@ -381,23 +389,63 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   /* Compute area */
   /* eqn. (7) */
   Anorm = 0.0f;
-  for (k = 0; k < 3; k++) {
-    /* we add a minus sign since dx is pi->x - pj->x */
-    A[k] = -Vi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) * wi *
-               hi_inv_dim -
-           Vj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) * wj *
-               hj_inv_dim;
-    Anorm += A[k] * A[k];
+  if (pi->density.wcorr > const_gizmo_min_wcorr &&
+      pj->density.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
+       different
+       smoothing lengths), then the expressions below are more stable. */
+    float Xi = Vi;
+    float Xj = Vj;
+#ifdef GIZMO_VOLUME_CORRECTION
+    if (fabsf(Vi - Vj) / fminf(Vi, Vj) > 1.5 * hydro_dimension) {
+      Xi = (Vi * hj + Vj * hi) / (hi + hj);
+      Xj = Xi;
+    }
+#endif
+    for (k = 0; k < 3; k++) {
+      /* we add a minus sign since dx is pi->x - pj->x */
+      A[k] = -Xi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) *
+                 wj * hj_inv_dim -
+             Xj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) *
+                 wi * hi_inv_dim;
+      Anorm += A[k] * A[k];
+    }
+  } else {
+    /* ill condition gradient matrix: revert to SPH face area */
+    Anorm = -(hidp1 * Vi * Vi * wi_dx + hjdp1 * Vj * Vj * wj_dx) * ri;
+    A[0] = -Anorm * dx[0];
+    A[1] = -Anorm * dx[1];
+    A[2] = -Anorm * dx[2];
+    Anorm *= Anorm * r2;
   }
 
-  if (!Anorm) {
+  if (Anorm == 0.) {
     /* if the interface has no area, nothing happens and we return */
     /* continuing results in dividing by zero and NaN's... */
     return;
   }
 
-  /* compute the normal vector of the interface */
   Anorm = sqrtf(Anorm);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* For stability reasons, we do require A and dx to have opposite
+     directions (basically meaning that the surface normal for the surface
+     always points from particle i to particle j, as it would in a real
+     moving-mesh code). If not, our scheme is no longer upwind and hence can
+     become unstable. */
+  float dA_dot_dx = A[0] * dx[0] + A[1] * dx[1] + A[2] * dx[2];
+  /* In GIZMO, Phil Hopkins reverts to an SPH integration scheme if this
+     happens. We curently just ignore this case and display a message. */
+  const float rdim = pow_dimension(r);
+  if (dA_dot_dx > 1.e-6 * rdim) {
+    message("Ill conditioned gradient matrix (%g %g %g %g %g)!", dA_dot_dx,
+            Anorm, Vi, Vj, r);
+  }
+#endif
+
+  /* compute the normal vector of the interface */
   for (k = 0; k < 3; k++) n_unit[k] = A[k] / Anorm;
 
   /* Compute interface position (relative to pi, since we don't need the actual
@@ -436,43 +484,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   /* we don't need to rotate, we can use the unit vector in the Riemann problem
    * itself (see GIZMO) */
 
-  if (Wi[0] < 0.0f || Wj[0] < 0.0f || Wi[4] < 0.0f || Wj[4] < 0.0f) {
-    printf("mindt: %g\n", mindt);
-    printf("WL: %g %g %g %g %g\n", pi->primitives.rho, pi->primitives.v[0],
-           pi->primitives.v[1], pi->primitives.v[2], pi->primitives.P);
-#ifdef USE_GRADIENTS
-    printf("dWL: %g %g %g %g %g\n", dWi[0], dWi[1], dWi[2], dWi[3], dWi[4]);
-#endif
-    printf("gradWL[0]: %g %g %g\n", pi->primitives.gradients.rho[0],
-           pi->primitives.gradients.rho[1], pi->primitives.gradients.rho[2]);
-    printf("gradWL[1]: %g %g %g\n", pi->primitives.gradients.v[0][0],
-           pi->primitives.gradients.v[0][1], pi->primitives.gradients.v[0][2]);
-    printf("gradWL[2]: %g %g %g\n", pi->primitives.gradients.v[1][0],
-           pi->primitives.gradients.v[1][1], pi->primitives.gradients.v[1][2]);
-    printf("gradWL[3]: %g %g %g\n", pi->primitives.gradients.v[2][0],
-           pi->primitives.gradients.v[2][1], pi->primitives.gradients.v[2][2]);
-    printf("gradWL[4]: %g %g %g\n", pi->primitives.gradients.P[0],
-           pi->primitives.gradients.P[1], pi->primitives.gradients.P[2]);
-    printf("WL': %g %g %g %g %g\n", Wi[0], Wi[1], Wi[2], Wi[3], Wi[4]);
-    printf("WR: %g %g %g %g %g\n", pj->primitives.rho, pj->primitives.v[0],
-           pj->primitives.v[1], pj->primitives.v[2], pj->primitives.P);
-#ifdef USE_GRADIENTS
-    printf("dWR: %g %g %g %g %g\n", dWj[0], dWj[1], dWj[2], dWj[3], dWj[4]);
-#endif
-    printf("gradWR[0]: %g %g %g\n", pj->primitives.gradients.rho[0],
-           pj->primitives.gradients.rho[1], pj->primitives.gradients.rho[2]);
-    printf("gradWR[1]: %g %g %g\n", pj->primitives.gradients.v[0][0],
-           pj->primitives.gradients.v[0][1], pj->primitives.gradients.v[0][2]);
-    printf("gradWR[2]: %g %g %g\n", pj->primitives.gradients.v[1][0],
-           pj->primitives.gradients.v[1][1], pj->primitives.gradients.v[1][2]);
-    printf("gradWR[3]: %g %g %g\n", pj->primitives.gradients.v[2][0],
-           pj->primitives.gradients.v[2][1], pj->primitives.gradients.v[2][2]);
-    printf("gradWR[4]: %g %g %g\n", pj->primitives.gradients.P[0],
-           pj->primitives.gradients.P[1], pj->primitives.gradients.P[2]);
-    printf("WR': %g %g %g %g %g\n", Wj[0], Wj[1], Wj[2], Wj[3], Wj[4]);
-    error("Negative density or pressure!\n");
-  }
-
   float totflux[5];
   riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);
 

src/hydro/Gizmo/hydro_io.h

@@ -127,7 +127,7 @@ float convert_Etot(struct engine* e, struct part* p) {
 void hydro_write_particles(struct part* parts, struct io_props* list,
                            int* num_fields) {
 
-  *num_fields = 14;
+  *num_fields = 11;
 
   /* List what we want to write */
   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -143,22 +143,16 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
                                               parts, primitives.P, convert_u);
   list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
                                  UNIT_CONV_NO_UNITS, parts, id);
-  list[6] = io_make_output_field("Acceleration", FLOAT, 3,
-                                 UNIT_CONV_ACCELERATION, parts, a_hydro);
-  list[7] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
+  list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
                                  primitives.rho);
-  list[8] = io_make_output_field("Volume", FLOAT, 1, UNIT_CONV_VOLUME, parts,
-                                 geometry.volume);
-  list[9] = io_make_output_field("GradDensity", FLOAT, 3, UNIT_CONV_DENSITY,
-                                 parts, primitives.gradients.rho);
-  list[10] = io_make_output_field_convert_part(
+  list[7] = io_make_output_field_convert_part(
       "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, primitives.P, convert_A);
-  list[11] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
-                                  parts, primitives.P);
-  list[12] =
+  list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
+                                 parts, primitives.P);
+  list[9] =
       io_make_output_field_convert_part("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
                                         parts, conserved.energy, convert_Etot);
-  list[13] = io_make_output_field("GravAcceleration", FLOAT, 3,
+  list[10] = io_make_output_field("GravAcceleration", FLOAT, 3,
                                   UNIT_CONV_ACCELERATION, parts, gravity.old_a);
 }
 

src/hydro/Gizmo/hydro_part.h

@@ -148,6 +148,9 @@ struct part {
     /* Total surface area of the particle. */
     float Atot;
 
+    /* Centroid of the "cell". */
+    float centroid[3];
+
   } geometry;
 
   /* Variables used for timestep calculation (currently not used). */
@@ -201,6 +204,8 @@ struct part {
     /* Previous value of the gravitational acceleration. */
     float old_a[3];
 
+    float grad_a[3][3];
+
     /* Previous value of the mass flux vector. */
     float old_mflux[3];
 

src/hydro/Gizmo/hydro_slope_limiters_face.h

@@ -53,14 +53,22 @@ hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
   if ((phimax + delta1) * phimax > 0.0f) {
     phiplus = phimax + delta1;
   } else {
-    phiplus = phimax / (1.0f + delta1 / fabs(phimax));
+    if (phimax != 0.) {
+      phiplus = phimax / (1.0f + delta1 / fabs(phimax));
+    } else {
+      phiplus = 0.;
+    }
   }
 
   /* if sign(phimin-delta1) == sign(phimin) */
   if ((phimin - delta1) * phimin > 0.0f) {
     phiminus = phimin - delta1;
   } else {
-    phiminus = phimin / (1.0f + delta1 / fabs(phimin));
+    if (phimin != 0.) {
+      phiminus = phimin / (1.0f + delta1 / fabs(phimin));
+    } else {
+      phiminus = 0.;
+    }
   }
 
   if (phi_i < phi_j) {

src/hydro/Gizmo/hydro_unphysical.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (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_UNPHYSICAL_H
+#define SWIFT_HYDRO_UNPHYSICAL_H
+
+#if defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
+
+#if defined(GIZMO_UNPHYSICAL_ERROR)
+
+/*! @brief Crash whenever an unphysical value is detected. */
+#define gizmo_unphysical_message(name, quantity) \
+  error("Unphysical " name " detected (%g)!", quantity);
+
+#elif defined(GIZMO_UNPHYSICAL_WARNING)
+
+/*! @brief Show a warning whenever an unphysical value is detected. */
+#define gizmo_unphysical_message(name, quantity) \
+  message("Unphysical " name " detected (%g), reset to 0!", quantity);
+
+#else
+
+/*! @brief Don't tell anyone an unphysical value was detected. */
+#define gizmo_unphysical_message(name, quantity)
+
+#endif
+
+#define gizmo_check_physical_quantity(name, quantity) \
+  if (quantity < 0.) {                                \
+    gizmo_unphysical_message(name, quantity);         \
+    quantity = 0.;                                    \
+  }
+
+#else  // defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
+
+#define gizmo_check_physical_quantity(name, quantity)
+
+#endif  // defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
+
+#endif  // SWIFT_HYDRO_UNPHYSICAL_H

src/hydro/Gizmo/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.;
+  p->v[1] = 0.;
+  p->v[2] = 0.;
+#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 particel data to act upon.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_velocities_prepare_force(struct part* restrict p,
+                               const struct xpart* restrict xp) {
+
+#ifndef GIZMO_FIX_PARTICLES
+  p->force.v_full[0] = xp->v_full[0];
+  p->force.v_full[1] = xp->v_full[1];
+  p->force.v_full[2] = xp->v_full[2];
+#endif
+}
+
+/**
+ * @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.;
+  p->v[1] = 0.;
+  p->v[2] = 0.;
+
+#else  // GIZMO_FIX_PARTICLES
+
+  if (p->conserved.mass > 0. && p->primitives.rho > 0.) {
+    /* Normal case: set particle velocity to fluid velocity. */
+    p->v[0] = p->conserved.momentum[0] / p->conserved.mass;
+    p->v[1] = p->conserved.momentum[1] / p->conserved.mass;
+    p->v[2] = p->conserved.momentum[2] / p->conserved.mass;
+  } else {
+    /* Vacuum particles have no fluid velocity. */
+    p->v[0] = 0.;
+    p->v[1] = 0.;
+    p->v[2] = 0.;
+  }
+
+#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.25;
+  const float etaR = eta * R;
+  const float xi = 1.;
+  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.9 * etaR) {
+    float fac = xi * soundspeed / d;
+    if (d < 1.1 * etaR) {
+      fac *= 5. * (d - 0.9 * etaR) / etaR;
+    }
+    p->v[0] -= ds[0] * fac;
+    p->v[1] -= ds[1] * fac;
+    p->v[2] -= ds[2] * fac;
+  }
+
+#endif  // GIZMO_STEER_MOTION
+
+#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.h

@@ -25,10 +25,8 @@
 #if defined(RIEMANN_SOLVER_EXACT)
 
 #define RIEMANN_SOLVER_IMPLEMENTATION "Exact Riemann solver (Toro 2009)"
-#if defined(EOS_IDEAL_GAS)
+#if defined(EOS_IDEAL_GAS) || defined(EOS_ISOTHERMAL_GAS)
 #include "riemann/riemann_exact.h"
-#elif defined(EOS_ISOTHERMAL_GAS)
-#include "riemann/riemann_exact_isothermal.h"
 #else
 #error "The Exact Riemann solver is incompatible with this equation of state!"
 #endif