diff --git a/examples/DiscPatch/HydroStatic/plot.py b/examples/DiscPatch/HydroStatic/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..2de749f9e3b3c287390218e09ea347d660f9ce8a
--- /dev/null
+++ b/examples/DiscPatch/HydroStatic/plot.py
@@ -0,0 +1,103 @@
+################################################################################
+# 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/>.
+#
+################################################################################
+
+##
+# This script plots the Disc-Patch_*.hdf5 snapshots.
+# It takes two (optional) parameters: the counter value of the first and last
+# snapshot to plot (default: 0 81).
+##
+
+import numpy as np
+import h5py
+import matplotlib
+matplotlib.use("Agg")
+import pylab as pl
+import glob
+import sys
+
+# Parameters
+surface_density = 10.
+scale_height = 100.
+z_disc = 200.
+utherm = 20.2615290634
+gamma = 5. / 3.
+
+start = 0
+stop = 81
+if len(sys.argv) > 1:
+  start = int(sys.argv[1])
+if len(sys.argv) > 2:
+  stop = int(sys.argv[2])
+
+# Get the analytic solution for the density
+def get_analytic_density(x):
+  return 0.5 * surface_density / scale_height / \
+           np.cosh( (x - z_disc) / scale_height )**2
+
+# Get the analytic solution for the (isothermal) pressure
+def get_analytic_pressure(x):
+  return (gamma - 1.) * utherm * get_analytic_density(x)
+
+# Get the data fields to plot from the snapshot file with the given name:
+#  snapshot time, z-coord, density, pressure, velocity norm
+def get_data(name):
+  file = h5py.File(name, "r")
+  coords = np.array(file["/PartType0/Coordinates"])
+  rho = np.array(file["/PartType0/Density"])
+  u = np.array(file["/PartType0/InternalEnergy"])
+  v = np.array(file["/PartType0/Velocities"])
+
+  P = (gamma - 1.) * rho * u
+
+  vtot = np.sqrt( v[:,0]**2 + v[:,1]**2 + v[:,2]**2 )
+
+  return float(file["/Header"].attrs["Time"]), coords[:,2], rho, P, vtot
+
+# scan the folder for snapshot files and plot all of them (within the requested
+# range)
+for f in sorted(glob.glob("Disc-Patch_*.hdf5")):
+  num = int(f[-8:-5])
+  if num < start or num > stop:
+    continue
+
+  print "processing", f, "..."
+
+  zrange = np.linspace(0., 400., 1000)
+  time, z, rho, P, v = get_data(f)
+
+  fig, ax = pl.subplots(3, 1, sharex = True)
+
+  ax[0].plot(z, rho, "r.")
+  ax[0].plot(zrange, get_analytic_density(zrange), "k-")
+  ax[0].set_ylabel("density")
+
+  ax[1].plot(z, v, "r.")
+  ax[1].plot(zrange, np.zeros(len(zrange)), "k-")
+  ax[1].set_ylabel("velocity norm")
+
+  ax[2].plot(z, P, "r.")
+  ax[2].plot(zrange, get_analytic_pressure(zrange), "k-")
+  ax[2].set_xlim(0., 400.)
+  ax[2].set_xlabel("z")
+  ax[2].set_ylabel("pressure")
+
+  pl.suptitle("t = {0:.2f}".format(time))
+
+  pl.savefig("{name}.png".format(name = f[:-5]))
+  pl.close()
diff --git a/examples/analyse_tasks.py b/examples/analyse_tasks.py
new file mode 100755
index 0000000000000000000000000000000000000000..04cd59feedba7ee41621ac0891d544c4aa294543
--- /dev/null
+++ b/examples/analyse_tasks.py
@@ -0,0 +1,183 @@
+#!/usr/bin/env python
+"""
+Usage:
+    analsyse_tasks.py [options] input.dat
+
+where input.dat is a thread info file for a step.  Use the '-y interval' flag
+of the swift command to create these.
+
+The output is an analysis of the task timings, including deadtime per thread
+and step, total amount of time spent for each task type, for the whole step
+and per thread and the minimum and maximum times spent per task type.
+
+This file is part of SWIFT.
+Copyright (c) 2017 Peter W. Draper (p.w.draper@durham.ac.uk)
+
+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/>.
+"""
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.collections as collections
+import matplotlib.ticker as plticker
+import pylab as pl
+import sys
+import argparse
+
+#  Handle the command line.
+parser = argparse.ArgumentParser(description="Analyse task dumps")
+
+parser.add_argument("input", help="Thread data file (-y output)")
+parser.add_argument("-v", "--verbose", dest="verbose",
+                    help="Verbose output (default: False)",
+                    default=False, action="store_true")
+
+args = parser.parse_args()
+infile = args.input
+
+#  Tasks and subtypes. Indexed as in tasks.h.
+TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
+             "init_grav", "ghost", "extra_ghost", "drift_part",
+             "drift_gpart", "kick1", "kick2", "timestep", "send", "recv",
+             "grav_top_level", "grav_long_range", "grav_mm", "grav_down",
+             "cooling", "sourceterms", "count"]
+
+SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
+            "tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
+
+#  Read input.
+data = pl.loadtxt( infile )
+
+maxthread = int(max(data[:,0])) + 1
+print "# Maximum thread id:", maxthread
+
+#  Recover the start and end time
+full_step = data[0,:]
+tic_step = int(full_step[4])
+toc_step = int(full_step[5])
+CPU_CLOCK = float(full_step[-1]) / 1000.0
+data = data[1:,:]
+if args.verbose:
+    print "CPU frequency:", CPU_CLOCK * 1000.0
+
+#  Avoid start and end times of zero.
+data = data[data[:,4] != 0]
+data = data[data[:,5] != 0]
+
+#  Calculate the time range.
+total_t = (toc_step - tic_step)/ CPU_CLOCK
+print "# Data range: ", total_t, "ms"
+
+#  Correct times to relative values.
+start_t = float(tic_step)
+data[:,4] -= start_t
+data[:,5] -= start_t
+
+tasks = {}
+tasks[-1] = []
+for i in range(maxthread):
+    tasks[i] = []
+
+#  Gather into by thread data.
+num_lines = pl.size(data) / 10
+for line in range(num_lines):
+    thread = int(data[line,0])
+    tic = int(data[line,4]) / CPU_CLOCK
+    toc = int(data[line,5]) / CPU_CLOCK
+    tasktype = int(data[line,1])
+    subtype = int(data[line,2])
+
+    tasks[thread].append([tic,toc,tasktype,subtype])
+
+#  Sort by tic and gather used thread ids.
+threadids = []
+for i in range(maxthread):
+    if len(tasks[i]) > 0:
+        tasks[i] = sorted(tasks[i], key=lambda task: task[0])
+        threadids.append(i)
+
+#  Times per task.
+print "# Task times:"
+print "# {0:<16s}: {1:>7s} {2:>9s} {3:>9s} {4:>9s} {5:>9s} {6:>9s}"\
+      .format("type/subtype", "count","minimum", "maximum",
+              "sum", "mean", "percent")
+alltasktimes = {}
+for i in threadids:
+    tasktimes = {}
+    for task in tasks[i]:
+        key = TASKTYPES[task[2]] + "/" + SUBTYPES[task[3]]
+        dt = task[1] - task[0]
+        if not key in tasktimes:
+            tasktimes[key] = []
+        tasktimes[key].append(dt)
+
+        if not key in alltasktimes:
+            alltasktimes[key] = []
+        alltasktimes[key].append(dt)
+
+    print "# Thread : ", i
+    for key in sorted(tasktimes.keys()):
+        taskmin = min(tasktimes[key])
+        taskmax = max(tasktimes[key])
+        tasksum = sum(tasktimes[key])
+        print "{0:18s}: {1:7d} {2:9.4f} {3:9.4f} {4:9.4f} {5:9.4f} {6:9.2f}"\
+              .format(key, len(tasktimes[key]), taskmin, taskmax, tasksum,
+                      tasksum / len(tasktimes[key]), tasksum / total_t * 100.0)
+    print
+
+print "# All threads : "
+for key in sorted(alltasktimes.keys()):
+    taskmin = min(alltasktimes[key])
+    taskmax = max(alltasktimes[key])
+    tasksum = sum(alltasktimes[key])
+    print "{0:18s}: {1:7d} {2:9.4f} {3:9.4f} {4:9.4f} {5:9.4f} {6:9.2f}"\
+          .format(key, len(alltasktimes[key]), taskmin, taskmax, tasksum,
+                  tasksum / len(alltasktimes[key]),
+                  tasksum / (len(threadids) * total_t) * 100.0)
+print
+
+#  Dead times.
+print "# Deadtimes:"
+print "# no.    : {0:>9s} {1:>9s} {2:>9s} {3:>9s} {4:>9s} {5:>9s}"\
+      .format("count", "minimum", "maximum", "sum", "mean", "percent")
+alldeadtimes = []
+for i in threadids:
+    deadtimes = []
+    last = 0
+    for task in tasks[i]:
+        dt = task[0] - last
+        deadtimes.append(dt)
+        last = task[1]
+    dt = total_t - last
+    deadtimes.append(dt)
+
+    deadmin = min(deadtimes)
+    deadmax = max(deadtimes)
+    deadsum = sum(deadtimes)
+    print "thread {0:2d}: {1:9d} {2:9.4f} {3:9.4f} {4:9.4f} {5:9.4f} {6:9.2f}"\
+          .format(i, len(deadtimes), deadmin, deadmax, deadsum,
+                  deadsum / len(deadtimes), deadsum / total_t * 100.0)
+    alldeadtimes.extend(deadtimes)
+
+deadmin = min(alldeadtimes)
+deadmax = max(alldeadtimes)
+deadsum = sum(alldeadtimes)
+print "all      : {0:9d} {1:9.4f} {2:9.4f} {3:9.4f} {4:9.4f} {5:9.2f}"\
+      .format(len(alldeadtimes), deadmin, deadmax, deadsum,
+              deadsum / len(alldeadtimes),
+              deadsum / (len(threadids) * total_t ) * 100.0)
+print
+
+
+sys.exit(0)
diff --git a/examples/analyse_tasks_MPI.py b/examples/analyse_tasks_MPI.py
new file mode 100755
index 0000000000000000000000000000000000000000..9feffaf67ec393257d75428e310a2e8b807df39a
--- /dev/null
+++ b/examples/analyse_tasks_MPI.py
@@ -0,0 +1,197 @@
+#!/usr/bin/env python
+"""
+Usage:
+    analsyse_tasks_MPI.py [options] input.dat
+
+where input.dat is a thread info file for an MPI step.  Use the '-y interval'
+flag of the swift command to create these.
+
+The output is an analysis of the task timings, including deadtime per thread
+and step, total amount of time spent for each task type, for the whole step
+and per thread and the minimum and maximum times spent per task type.
+
+This file is part of SWIFT.
+Copyright (c) 2017 Peter W. Draper (p.w.draper@durham.ac.uk)
+
+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/>.
+"""
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.collections as collections
+import matplotlib.ticker as plticker
+import pylab as pl
+import sys
+import argparse
+
+#  Handle the command line.
+parser = argparse.ArgumentParser(description="Analyse task dumps")
+
+parser.add_argument("input", help="Thread data file (-y output)")
+parser.add_argument("-v", "--verbose", dest="verbose",
+                    help="Verbose output (default: False)",
+                    default=False, action="store_true")
+
+args = parser.parse_args()
+infile = args.input
+
+#  Tasks and subtypes. Indexed as in tasks.h.
+TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
+             "init_grav", "ghost", "extra_ghost", "drift_part",
+             "drift_gpart", "kick1", "kick2", "timestep", "send", "recv",
+             "grav_top_level", "grav_long_range", "grav_mm", "grav_down",
+             "cooling", "sourceterms", "count"]
+
+SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
+            "tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
+
+#  Read input.
+data = pl.loadtxt( infile )
+
+#  Get the CPU clock to convert ticks into milliseconds.
+full_step = data[0,:]
+CPU_CLOCK = float(full_step[-1]) / 1000.0
+if args.verbose:
+    print "# CPU frequency:", CPU_CLOCK * 1000.0
+
+nranks = int(max(data[:,0])) + 1
+print "# Number of ranks:", nranks
+maxthread = int(max(data[:,1])) + 1
+print "# Maximum thread id:", maxthread
+
+#  Avoid start and end times of zero.
+sdata = data[data[:,5] != 0]
+sdata = data[data[:,6] != 0]
+
+#  Now we process all the ranks.
+for rank in range(nranks):
+    print "# Rank", rank
+    data = sdata[sdata[:,0] == rank]
+
+    #  Recover the start and end time
+    full_step = data[0,:]
+    tic_step = int(full_step[5])
+    toc_step = int(full_step[6])
+    data = data[1:,:]
+
+    #  Avoid start and end times of zero.
+    data = data[data[:,5] != 0]
+    data = data[data[:,6] != 0]
+
+    #  Calculate the time range.
+    total_t = (toc_step - tic_step)/ CPU_CLOCK
+    print "# Data range: ", total_t, "ms"
+
+    #  Correct times to relative values.
+    start_t = float(tic_step)
+    data[:,5] -= start_t
+    data[:,6] -= start_t
+    end_t = (toc_step - start_t) / CPU_CLOCK
+
+    tasks = {}
+    tasks[-1] = []
+    for i in range(maxthread):
+        tasks[i] = []
+
+    #  Gather into by thread data.
+    num_lines = pl.size(data) / 12
+    for line in range(num_lines):
+        thread = int(data[line,1])
+        tic = int(data[line,5]) / CPU_CLOCK
+        toc = int(data[line,6]) / CPU_CLOCK
+        tasktype = int(data[line,2])
+        subtype = int(data[line,3])
+
+        tasks[thread].append([tic,toc,tasktype,subtype])
+
+    #  Sort by tic and gather used threads.
+    threadids = []
+    for i in range(maxthread):
+        tasks[i] = sorted(tasks[i], key=lambda task: task[0])
+        threadids.append(i)
+
+    #  Times per task.
+    print "# Task times:"
+    print "# {0:<16s}: {1:>7s} {2:>9s} {3:>9s} {4:>9s} {5:>9s} {6:>9s}"\
+          .format("type/subtype", "count","minimum", "maximum",
+                  "sum", "mean", "percent")
+    alltasktimes = {}
+    for i in threadids:
+        tasktimes = {}
+        for task in tasks[i]:
+            key = TASKTYPES[task[2]] + "/" + SUBTYPES[task[3]]
+            dt = task[1] - task[0]
+            if not key in tasktimes:
+                tasktimes[key] = []
+            tasktimes[key].append(dt)
+
+            if not key in alltasktimes:
+                alltasktimes[key] = []
+            alltasktimes[key].append(dt)
+
+        print "# Thread : ", i
+        for key in sorted(tasktimes.keys()):
+            taskmin = min(tasktimes[key])
+            taskmax = max(tasktimes[key])
+            tasksum = sum(tasktimes[key])
+            print "{0:18s}: {1:7d} {2:9.4f} {3:9.4f} {4:9.4f} {5:9.4f} {6:9.2f}"\
+                  .format(key, len(tasktimes[key]), taskmin, taskmax, tasksum,
+                          tasksum / len(tasktimes[key]), tasksum / total_t * 100.0)
+        print
+
+    print "# All threads : "
+    for key in sorted(alltasktimes.keys()):
+        taskmin = min(alltasktimes[key])
+        taskmax = max(alltasktimes[key])
+        tasksum = sum(alltasktimes[key])
+        print "{0:18s}: {1:7d} {2:9.4f} {3:9.4f} {4:9.4f} {5:9.4f} {6:9.2f}"\
+              .format(key, len(alltasktimes[key]), taskmin, taskmax, tasksum,
+                      tasksum / len(alltasktimes[key]),
+                      tasksum / (len(threadids) * total_t) * 100.0)
+    print
+
+    #  Dead times.
+    print "# Deadtimes:"
+    print "# no.    : {0:>9s} {1:>9s} {2:>9s} {3:>9s} {4:>9s} {5:>9s}"\
+          .format("count", "minimum", "maximum", "sum", "mean", "percent")
+    alldeadtimes = []
+    for i in threadids:
+        deadtimes = []
+        last = 0
+        for task in tasks[i]:
+            dt = task[0] - last
+            deadtimes.append(dt)
+            last = task[1]
+        dt = total_t - last
+        deadtimes.append(dt)
+
+        deadmin = min(deadtimes)
+        deadmax = max(deadtimes)
+        deadsum = sum(deadtimes)
+        print "thread {0:2d}: {1:9d} {2:9.4f} {3:9.4f} {4:9.4f} {5:9.4f} {6:9.2f}"\
+              .format(i, len(deadtimes), deadmin, deadmax, deadsum,
+                      deadsum / len(deadtimes), deadsum / total_t * 100.0)
+        alldeadtimes.extend(deadtimes)
+
+    deadmin = min(alldeadtimes)
+    deadmax = max(alldeadtimes)
+    deadsum = sum(alldeadtimes)
+    print "all      : {0:9d} {1:9.4f} {2:9.4f} {3:9.4f} {4:9.4f} {5:9.2f}"\
+          .format(len(alldeadtimes), deadmin, deadmax, deadsum,
+                  deadsum / len(alldeadtimes),
+              deadsum / (len(threadids) * total_t ) * 100.0)
+    print
+
+
+sys.exit(0)
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 14bc60bc1e1c05ecdc66fb7ac828102b1d5748bf..8006c1a325845d6e9fec655b809310a63daa9ddb 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -107,6 +107,12 @@ DiscPatchPotential:
   timestep_mult:   0.03     # Dimensionless pre-factor for the time-step condition
   growth_time:     5.       # (Optional) Time for the disc to grow to its final size (multiple of the dynamical time)
 
+# Sine Wave potential
+SineWavePotential:
+  amplitude:        10.     # Amplitude of the sine wave (internal units)
+  timestep_limit:   1.      # Time-step dimensionless pre-factor.
+  growth_time:      0.      # (Optional) Time for the potential to grow to its final size.
+ 
 # Parameters related to cooling function  ----------------------------------------------
 
 # Constant du/dt cooling function
diff --git a/examples/plot_tasks.py b/examples/plot_tasks.py
index 75ac2db65cac748606c143b3c691a167b314a926..88f176687db8116cfd4370970769164985e4d366 100755
--- a/examples/plot_tasks.py
+++ b/examples/plot_tasks.py
@@ -89,9 +89,10 @@ pl.rcParams.update(PLOT_PARAMS)
 
 #  Tasks and subtypes. Indexed as in tasks.h.
 TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
-             "init_grav", "ghost", "extra_ghost", "drift", "kick1", "kick2",
-             "timestep", "send", "recv", "grav_top_level", "grav_long_range",
-             "grav_mm", "grav_down", "cooling", "sourceterms", "count"]
+             "init_grav", "ghost", "extra_ghost", "drift_part",
+             "drift_gpart", "kick1", "kick2", "timestep", "send", "recv",
+             "grav_top_level", "grav_long_range", "grav_mm", "grav_down",
+             "cooling", "sourceterms", "count"]
 
 SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
             "tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
@@ -105,14 +106,14 @@ FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
 
 #  A number of colours for the various types. Recycled when there are
 #  more task types than colours...
-
 colours = ["cyan", "lightgray", "darkblue", "yellow", "tan", "dodgerblue",
-           "sienna", "aquamarine", "bisque", "blue", "green", "brown",
-           "purple", "mocassin", "olivedrab", "chartreuse", "darksage",
-           "darkgreen", "green", "mediumseagreen", "mediumaquamarine",
-           "darkslategrey", "mediumturquoise", "black", "cadetblue", "skyblue",
-           "red", "slategray", "gold", "slateblue", "blueviolet",
-           "mediumorchid", "firebrick", "magenta", "hotpink", "pink"]
+           "sienna", "aquamarine", "bisque", "blue", "green", "lightgreen",
+           "brown", "purple", "moccasin", "olivedrab", "chartreuse",
+           "darksage", "darkgreen", "green", "mediumseagreen",
+           "mediumaquamarine", "darkslategrey", "mediumturquoise",
+           "black", "cadetblue", "skyblue", "red", "slategray", "gold",
+           "slateblue", "blueviolet", "mediumorchid", "firebrick",
+           "magenta", "hotpink", "pink", "orange", "lightgreen"]
 maxcolours = len(colours)
 
 #  Set colours of task/subtype.
@@ -134,9 +135,9 @@ for task in SUBTYPES:
 #  For fiddling with colours...
 if args.verbose:
     print "#Selected colours:"
-    for task in TASKCOLOURS.keys():
+    for task in sorted(TASKCOLOURS.keys()):
         print "# " + task + ": " + TASKCOLOURS[task]
-    for task in SUBCOLOURS.keys():
+    for task in sorted(SUBCOLOURS.keys()):
         print "# " + task + ": " + SUBCOLOURS[task]
 
 #  Read input.
@@ -161,7 +162,7 @@ data = data[data[:,5] != 0]
 #  Calculate the time range, if not given.
 delta_t = delta_t * CPU_CLOCK
 if delta_t == 0:
-    dt = max(data[:,5]) - min(data[:,4])
+    dt = toc_step - tic_step
     if dt > delta_t:
         delta_t = dt
     print "Data range: ", delta_t / CPU_CLOCK, "ms"
diff --git a/examples/plot_tasks_MPI.py b/examples/plot_tasks_MPI.py
index a1f81bb1789c73ebe21391509c7d7674f6e43100..83465aee87e8b641775d760fa4db2f06b125dd8b 100755
--- a/examples/plot_tasks_MPI.py
+++ b/examples/plot_tasks_MPI.py
@@ -95,9 +95,10 @@ pl.rcParams.update(PLOT_PARAMS)
 
 #  Tasks and subtypes. Indexed as in tasks.h.
 TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
-             "init_grav", "ghost", "extra_ghost", "drift", "kick1", "kick2",
-             "timestep", "send", "recv", "grav_top_level", "grav_long_range",
-             "grav_mm", "grav_down", "cooling", "sourceterms", "count"]
+             "init_grav", "ghost", "extra_ghost", "drift_part", "drift_gpart",
+             "kick1", "kick2", "timestep", "send", "recv", "grav_top_level",
+             "grav_long_range", "grav_mm", "grav_down", "cooling",
+             "sourceterms", "count"]
 
 SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
             "tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
@@ -111,15 +112,14 @@ FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
 
 #  A number of colours for the various types. Recycled when there are
 #  more task types than colours...
-
 colours = ["cyan", "lightgray", "darkblue", "yellow", "tan", "dodgerblue",
-           "sienna", "aquamarine", "bisque", "blue", "green", "brown",
-           "purple", "mocassin", "olivedrab", "chartreuse", "darksage",
-           "darkgreen", "green", "mediumseagreen", "mediumaquamarine",
-           "darkslategrey", "mediumturquoise", "black", "cadetblue", "skyblue",
-           "red", "slategray", "gold", "slateblue", "blueviolet",
-           "mediumorchid", "firebrick", "magenta", "hotpink", "pink",
-           "orange", "lightgreen"]
+           "sienna", "aquamarine", "bisque", "blue", "green", "lightgreen",
+           "brown", "purple", "moccasin", "olivedrab", "chartreuse",
+           "darksage", "darkgreen", "green", "mediumseagreen",
+           "mediumaquamarine", "darkslategrey", "mediumturquoise",
+           "black", "cadetblue", "skyblue", "red", "slategray", "gold",
+           "slateblue", "blueviolet", "mediumorchid", "firebrick",
+           "magenta", "hotpink", "pink", "orange", "lightgreen"]
 maxcolours = len(colours)
 
 #  Set colours of task/subtype.
@@ -141,9 +141,9 @@ for task in SUBTYPES:
 #  For fiddling with colours...
 if args.verbose:
     print "#Selected colours:"
-    for task in TASKCOLOURS.keys():
+    for task in sorted(TASKCOLOURS.keys()):
         print "# " + task + ": " + TASKCOLOURS[task]
-    for task in SUBCOLOURS.keys():
+    for task in sorted(SUBCOLOURS.keys()):
         print "# " + task + ": " + SUBCOLOURS[task]
 
 #  Read input.
@@ -171,12 +171,14 @@ delta_t = delta_t * CPU_CLOCK
 if delta_t == 0:
     for rank in range(nranks):
         data = sdata[sdata[:,0] == rank]
-        dt = max(data[:,6]) - min(data[:,5])
+        full_step = data[0,:]
+        tic_step = int(full_step[5])
+        toc_step = int(full_step[6])
+        dt = toc_step - tic_step
         if dt > delta_t:
             delta_t = dt
     print "Data range: ", delta_t / CPU_CLOCK, "ms"
 
-
 # Once more doing the real gather and plots this time.
 for rank in range(nranks):
     data = sdata[sdata[:,0] == rank]
@@ -186,6 +188,8 @@ for rank in range(nranks):
     tic_step = int(full_step[5])
     toc_step = int(full_step[6])
     data = data[1:,:]
+    typesseen = []
+    nethread = 0
 
     #  Dummy image for ranks that have no tasks.
     if data.size == 0:
diff --git a/examples/process_plot_tasks b/examples/process_plot_tasks
index c33531269617e81d4e8f18f1528b43492b08cf26..b46fce03d8c5f21046a0e4a95a304e006c7b2293 100755
--- a/examples/process_plot_tasks
+++ b/examples/process_plot_tasks
@@ -57,6 +57,7 @@ done
 #  And process them,
 echo "Processing thread info files..."
 echo $list | xargs -P $NPROCS -n 3 /bin/bash -c "./plot_tasks.py --expand 1 --limit $TIMERANGE --width 16 --height 4 \$0 \$2 "
+echo $list | xargs -P $NPROCS -n 3 /bin/bash -c "./analyse_tasks.py \$0 > \$2.stats"
 
 echo "Writing output index.html file"
 #  Construct document - serial.
@@ -75,8 +76,21 @@ echo $list | xargs -n 3 | while read f s g; do
 <h2>Step $s</h2>
 EOF
     cat <<EOF >> index.html
-<a href="step${s}r${i}.png"><img src="step${s}r${i}.png" width=400px/></a>
+<a href="step${s}r${i}.html"><img src="step${s}r${i}.png" width=400px/></a>
 EOF
+    cat <<EOF > step${s}r${i}.html
+ <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<body>
+<img src="step${s}r${i}.png">
+<pre>
+EOF
+cat step${s}r${i}.stats >> step${s}r${i}.html
+cat <<EOF >> step${s}r${i}.html
+</body>
+</html>
+EOF
+
 done
 
 cat <<EOF >> index.html
diff --git a/examples/process_plot_tasks_MPI b/examples/process_plot_tasks_MPI
index 268a5393dc94f8b05cb4ff0a71c45b3b8e554c84..b2672b3711823eb87d0bede5b1ffd8945a735f98 100755
--- a/examples/process_plot_tasks_MPI
+++ b/examples/process_plot_tasks_MPI
@@ -62,6 +62,7 @@ nrank=$(($nrank-1))
 #  And process them,
 echo "Processing thread info files..."
 echo $list | xargs -P $NPROCS -n 3 /bin/bash -c "./plot_tasks_MPI.py --expand 1 --limit $TIMERANGE \$0 \$2 "
+echo $list | xargs -P $NPROCS -n 3 /bin/bash -c "./analyse_tasks_MPI.py \$0 > \$2.stats"
 
 echo "Writing output index.html file"
 #  Construct document - serial.
@@ -78,12 +79,31 @@ EOF
 echo $list | xargs -n 3 | while read f s g; do
     cat <<EOF >> index.html
 <h2>Step $s</h2>
+<ul style="list-style-type:none">
+<li>
 EOF
     for i in $(seq 0 $nrank); do
-        cat <<EOF >> index.html
-<a href="step${s}r${i}.png"><img src="step${s}r${i}.png" width=400px/></a>
-EOF
+        cat <<EOF2 >> index.html
+<a href="step${s}r${i}.html"><img src="step${s}r${i}.png" width=400px/></a>
+EOF2
+    cat <<EOF2 > step${s}r${i}.html
+ <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<body>
+<img src="step${s}r${i}.png">
+<pre>
+EOF2
+cat step${s}r.stats >> step${s}r${i}.html
+cat <<EOF2 >> step${s}r${i}.html
+</pre>
+</body>
+</html>
+EOF2
     done
+cat <<EOF >> index.html
+</li>
+</ul>
+EOF
 done
 
 cat <<EOF >> index.html
diff --git a/src/Makefile.am b/src/Makefile.am
index 809f5c34642392d001f1987c1c8926b8b97eae1e..2ddcdb0908201c65053d7cc5380a4217277b5c13 100644
--- a/src/Makefile.am
+++ b/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 \
diff --git a/src/cell.c b/src/cell.c
index 657cd85cea63635c532bd05e4425600e2f7fefb5..dbccfd2f42cabf38417cd87de0450489240884be 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1385,6 +1385,9 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
         scheduler_activate(s, ci->recv_xv);
         if (cell_is_active(ci, e)) {
           scheduler_activate(s, ci->recv_rho);
+#ifdef EXTRA_HYDRO_LOOP
+          scheduler_activate(s, ci->recv_gradient);
+#endif
           scheduler_activate(s, ci->recv_ti);
         }
 
@@ -1404,12 +1407,21 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
         if (t->type == task_type_pair) scheduler_activate(s, cj->drift_part);
 
         if (cell_is_active(cj, e)) {
+
           for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
                l = l->next)
             ;
           if (l == NULL) error("Missing link to send_rho task.");
           scheduler_activate(s, l->t);
 
+#ifdef EXTRA_HYDRO_LOOP
+          for (l = cj->send_gradient;
+               l != NULL && l->t->cj->nodeID != ci->nodeID; l = l->next)
+            ;
+          if (l == NULL) error("Missing link to send_gradient task.");
+          scheduler_activate(s, l->t);
+#endif
+
           for (l = cj->send_ti; l != NULL && l->t->cj->nodeID != ci->nodeID;
                l = l->next)
             ;
@@ -1423,6 +1435,9 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
         scheduler_activate(s, cj->recv_xv);
         if (cell_is_active(cj, e)) {
           scheduler_activate(s, cj->recv_rho);
+#ifdef EXTRA_HYDRO_LOOP
+          scheduler_activate(s, cj->recv_gradient);
+#endif
           scheduler_activate(s, cj->recv_ti);
         }
 
@@ -1442,12 +1457,21 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
         if (t->type == task_type_pair) scheduler_activate(s, ci->drift_part);
 
         if (cell_is_active(ci, e)) {
+
           for (l = ci->send_rho; l != NULL && l->t->cj->nodeID != cj->nodeID;
                l = l->next)
             ;
           if (l == NULL) error("Missing link to send_rho task.");
           scheduler_activate(s, l->t);
 
+#ifdef EXTRA_HYDRO_LOOP
+          for (l = ci->send_gradient;
+               l != NULL && l->t->cj->nodeID != cj->nodeID; l = l->next)
+            ;
+          if (l == NULL) error("Missing link to send_gradient task.");
+          scheduler_activate(s, l->t);
+#endif
+
           for (l = ci->send_ti; l != NULL && l->t->cj->nodeID != cj->nodeID;
                l = l->next)
             ;
diff --git a/src/const.h b/src/const.h
index 6962ee8bca32e92664e3f20cdb23e7cb6fbc4abd..141eb48acc633542aa98655caa8debdd2dbce530 100644
--- a/src/const.h
+++ b/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
diff --git a/src/engine.c b/src/engine.c
index 28b1d4b2e8b42cb865125f610f718b1d7130bd63..417c9f626d7e2f8d96d49d8d2bed942102b96e4f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2584,6 +2584,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         scheduler_activate(s, ci->recv_xv);
         if (cell_is_active(ci, e)) {
           scheduler_activate(s, ci->recv_rho);
+#ifdef EXTRA_HYDRO_LOOP
+          scheduler_activate(s, ci->recv_gradient);
+#endif
           scheduler_activate(s, ci->recv_ti);
         }
 
@@ -2603,12 +2606,21 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         if (t->type == task_type_pair) scheduler_activate(s, cj->drift_part);
 
         if (cell_is_active(cj, e)) {
+
           for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
                l = l->next)
             ;
           if (l == NULL) error("Missing link to send_rho task.");
           scheduler_activate(s, l->t);
 
+#ifdef EXTRA_HYDRO_LOOP
+          for (l = cj->send_gradient;
+               l != NULL && l->t->cj->nodeID != ci->nodeID; l = l->next)
+            ;
+          if (l == NULL) error("Missing link to send_gradient task.");
+          scheduler_activate(s, l->t);
+#endif
+
           for (l = cj->send_ti; l != NULL && l->t->cj->nodeID != ci->nodeID;
                l = l->next)
             ;
@@ -2622,6 +2634,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         scheduler_activate(s, cj->recv_xv);
         if (cell_is_active(cj, e)) {
           scheduler_activate(s, cj->recv_rho);
+#ifdef EXTRA_HYDRO_LOOP
+          scheduler_activate(s, cj->recv_gradient);
+#endif
           scheduler_activate(s, cj->recv_ti);
         }
 
@@ -2647,6 +2662,14 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
           if (l == NULL) error("Missing link to send_rho task.");
           scheduler_activate(s, l->t);
 
+#ifdef EXTRA_HYDRO_LOOP
+          for (l = ci->send_gradient;
+               l != NULL && l->t->cj->nodeID != cj->nodeID; l = l->next)
+            ;
+          if (l == NULL) error("Missing link to send_gradient task.");
+          scheduler_activate(s, l->t);
+#endif
+
           for (l = ci->send_ti; l != NULL && l->t->cj->nodeID != cj->nodeID;
                l = l->next)
             ;
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 2e340a03b99ae51bc49a2e57456f4d6838d62f21..6d39c54d2ddc3571ac34c54fc9eede6f7dee6ac5 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/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 */
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index a5c1e9038d0d3de6896afe773e3193a2304a6b6b..5ad6d87619a7629a703a8b9c03d089e69ffbdf7d 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/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
diff --git a/src/hydro/Gizmo/hydro_gradients_gizmo.h b/src/hydro/Gizmo/hydro_gradients_gizmo.h
index aa6e4406b94e7a5cafcd0ca556162476003477de..ee3ad6919f81f042ceacc5db8b4e818d63c90266 100644
--- a/src/hydro/Gizmo/hydro_gradients_gizmo.h
+++ b/src/hydro/Gizmo/hydro_gradients_gizmo.h
@@ -45,6 +45,18 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
   p->primitives.gradients.P[1] = 0.0f;
   p->primitives.gradients.P[2] = 0.0f;
 
+  p->gravity.grad_a[0][0] = 0.0f;
+  p->gravity.grad_a[0][1] = 0.0f;
+  p->gravity.grad_a[0][2] = 0.0f;
+
+  p->gravity.grad_a[1][0] = 0.0f;
+  p->gravity.grad_a[1][1] = 0.0f;
+  p->gravity.grad_a[1][2] = 0.0f;
+
+  p->gravity.grad_a[2][0] = 0.0f;
+  p->gravity.grad_a[2][1] = 0.0f;
+  p->gravity.grad_a[2][2] = 0.0f;
+
   hydro_slope_limit_cell_init(p);
 }
 
@@ -93,56 +105,146 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
   xi = r * hi_inv;
   kernel_deval(xi, &wi, &wi_dx);
 
-  /* Compute gradients for pi */
-  /* there is a sign difference w.r.t. eqn. (6) because of the inverse
-   * definition of dx */
-  pi->primitives.gradients.rho[0] +=
-      (Wi[0] - Wj[0]) * wi *
-      (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-  pi->primitives.gradients.rho[1] +=
-      (Wi[0] - Wj[0]) * wi *
-      (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-  pi->primitives.gradients.rho[2] +=
-      (Wi[0] - Wj[0]) * wi *
-      (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
-  pi->primitives.gradients.v[0][0] +=
-      (Wi[1] - Wj[1]) * wi *
-      (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-  pi->primitives.gradients.v[0][1] +=
-      (Wi[1] - Wj[1]) * wi *
-      (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-  pi->primitives.gradients.v[0][2] +=
-      (Wi[1] - Wj[1]) * wi *
-      (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-  pi->primitives.gradients.v[1][0] +=
-      (Wi[2] - Wj[2]) * wi *
-      (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-  pi->primitives.gradients.v[1][1] +=
-      (Wi[2] - Wj[2]) * wi *
-      (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-  pi->primitives.gradients.v[1][2] +=
-      (Wi[2] - Wj[2]) * wi *
-      (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-  pi->primitives.gradients.v[2][0] +=
-      (Wi[3] - Wj[3]) * wi *
-      (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-  pi->primitives.gradients.v[2][1] +=
-      (Wi[3] - Wj[3]) * wi *
-      (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-  pi->primitives.gradients.v[2][2] +=
-      (Wi[3] - Wj[3]) * wi *
-      (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
-  pi->primitives.gradients.P[0] +=
-      (Wi[4] - Wj[4]) * wi *
-      (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-  pi->primitives.gradients.P[1] +=
-      (Wi[4] - Wj[4]) * wi *
-      (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-  pi->primitives.gradients.P[2] +=
-      (Wi[4] - Wj[4]) * wi *
-      (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+  if (pi->density.wcorr > const_gizmo_min_wcorr) {
+    /* Compute gradients for pi */
+    /* there is a sign difference w.r.t. eqn. (6) because of the inverse
+     * definition of dx */
+    pi->primitives.gradients.rho[0] +=
+        (Wi[0] - Wj[0]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.rho[1] +=
+        (Wi[0] - Wj[0]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.rho[2] +=
+        (Wi[0] - Wj[0]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->primitives.gradients.v[0][0] +=
+        (Wi[1] - Wj[1]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.v[0][1] +=
+        (Wi[1] - Wj[1]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.v[0][2] +=
+        (Wi[1] - Wj[1]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+    pi->primitives.gradients.v[1][0] +=
+        (Wi[2] - Wj[2]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.v[1][1] +=
+        (Wi[2] - Wj[2]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.v[1][2] +=
+        (Wi[2] - Wj[2]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+    pi->primitives.gradients.v[2][0] +=
+        (Wi[3] - Wj[3]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.v[2][1] +=
+        (Wi[3] - Wj[3]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.v[2][2] +=
+        (Wi[3] - Wj[3]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->primitives.gradients.P[0] +=
+        (Wi[4] - Wj[4]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.P[1] +=
+        (Wi[4] - Wj[4]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.P[2] +=
+        (Wi[4] - Wj[4]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->gravity.grad_a[0][0] +=
+        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->gravity.grad_a[0][1] +=
+        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->gravity.grad_a[0][2] +=
+        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->gravity.grad_a[1][0] +=
+        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->gravity.grad_a[1][1] +=
+        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->gravity.grad_a[1][2] +=
+        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->gravity.grad_a[2][0] +=
+        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->gravity.grad_a[2][1] +=
+        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->gravity.grad_a[2][2] +=
+        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+  } else {
+    /* The gradient matrix was not well-behaved, switch to SPH gradients */
+
+    pi->primitives.gradients.rho[0] -=
+        wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+    pi->primitives.gradients.rho[1] -=
+        wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+    pi->primitives.gradients.rho[2] -=
+        wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+    pi->primitives.gradients.v[0][0] -=
+        wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+    pi->primitives.gradients.v[0][1] -=
+        wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+    pi->primitives.gradients.v[0][2] -=
+        wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+
+    pi->primitives.gradients.v[1][0] -=
+        wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+    pi->primitives.gradients.v[1][1] -=
+        wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+    pi->primitives.gradients.v[1][2] -=
+        wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+
+    pi->primitives.gradients.v[2][0] -=
+        wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+    pi->primitives.gradients.v[2][1] -=
+        wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+    pi->primitives.gradients.v[2][2] -=
+        wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+    pi->primitives.gradients.P[0] -=
+        wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+    pi->primitives.gradients.P[1] -=
+        wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+    pi->primitives.gradients.P[2] -=
+        wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+    pi->gravity.grad_a[0][0] -=
+        wi_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+    pi->gravity.grad_a[0][1] -=
+        wi_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+    pi->gravity.grad_a[0][2] -=
+        wi_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+
+    pi->gravity.grad_a[1][0] -=
+        wi_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+    pi->gravity.grad_a[1][1] -=
+        wi_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+    pi->gravity.grad_a[1][2] -=
+        wi_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+
+    pi->gravity.grad_a[2][0] -=
+        wi_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+    pi->gravity.grad_a[2][1] -=
+        wi_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+    pi->gravity.grad_a[2][2] -=
+        wi_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+  }
 
   hydro_slope_limit_cell_collect(pi, pj, r);
 
@@ -151,57 +253,146 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
   xj = r * hj_inv;
   kernel_deval(xj, &wj, &wj_dx);
 
-  /* Compute gradients for pj */
-  /* there is no sign difference w.r.t. eqn. (6) because dx is now what we
-   * want
-   * it to be */
-  pj->primitives.gradients.rho[0] +=
-      (Wi[0] - Wj[0]) * wj *
-      (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
-  pj->primitives.gradients.rho[1] +=
-      (Wi[0] - Wj[0]) * wj *
-      (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
-  pj->primitives.gradients.rho[2] +=
-      (Wi[0] - Wj[0]) * wj *
-      (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
-  pj->primitives.gradients.v[0][0] +=
-      (Wi[1] - Wj[1]) * wj *
-      (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
-  pj->primitives.gradients.v[0][1] +=
-      (Wi[1] - Wj[1]) * wj *
-      (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
-  pj->primitives.gradients.v[0][2] +=
-      (Wi[1] - Wj[1]) * wj *
-      (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-  pj->primitives.gradients.v[1][0] +=
-      (Wi[2] - Wj[2]) * wj *
-      (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
-  pj->primitives.gradients.v[1][1] +=
-      (Wi[2] - Wj[2]) * wj *
-      (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
-  pj->primitives.gradients.v[1][2] +=
-      (Wi[2] - Wj[2]) * wj *
-      (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-  pj->primitives.gradients.v[2][0] +=
-      (Wi[3] - Wj[3]) * wj *
-      (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
-  pj->primitives.gradients.v[2][1] +=
-      (Wi[3] - Wj[3]) * wj *
-      (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
-  pj->primitives.gradients.v[2][2] +=
-      (Wi[3] - Wj[3]) * wj *
-      (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
-  pj->primitives.gradients.P[0] +=
-      (Wi[4] - Wj[4]) * wj *
-      (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
-  pj->primitives.gradients.P[1] +=
-      (Wi[4] - Wj[4]) * wj *
-      (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
-  pj->primitives.gradients.P[2] +=
-      (Wi[4] - Wj[4]) * wj *
-      (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+  if (pj->density.wcorr > const_gizmo_min_wcorr) {
+    /* Compute gradients for pj */
+    /* there is no sign difference w.r.t. eqn. (6) because dx is now what we
+     * want
+     * it to be */
+    pj->primitives.gradients.rho[0] +=
+        (Wi[0] - Wj[0]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->primitives.gradients.rho[1] +=
+        (Wi[0] - Wj[0]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->primitives.gradients.rho[2] +=
+        (Wi[0] - Wj[0]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+    pj->primitives.gradients.v[0][0] +=
+        (Wi[1] - Wj[1]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->primitives.gradients.v[0][1] +=
+        (Wi[1] - Wj[1]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->primitives.gradients.v[0][2] +=
+        (Wi[1] - Wj[1]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+    pj->primitives.gradients.v[1][0] +=
+        (Wi[2] - Wj[2]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->primitives.gradients.v[1][1] +=
+        (Wi[2] - Wj[2]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->primitives.gradients.v[1][2] +=
+        (Wi[2] - Wj[2]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+    pj->primitives.gradients.v[2][0] +=
+        (Wi[3] - Wj[3]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->primitives.gradients.v[2][1] +=
+        (Wi[3] - Wj[3]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->primitives.gradients.v[2][2] +=
+        (Wi[3] - Wj[3]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+    pj->primitives.gradients.P[0] +=
+        (Wi[4] - Wj[4]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->primitives.gradients.P[1] +=
+        (Wi[4] - Wj[4]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->primitives.gradients.P[2] +=
+        (Wi[4] - Wj[4]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+    pj->gravity.grad_a[0][0] +=
+        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->gravity.grad_a[0][1] +=
+        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->gravity.grad_a[0][2] +=
+        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+    pj->gravity.grad_a[1][0] +=
+        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->gravity.grad_a[1][1] +=
+        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->gravity.grad_a[1][2] +=
+        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+    pj->gravity.grad_a[2][0] +=
+        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->gravity.grad_a[2][1] +=
+        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->gravity.grad_a[2][2] +=
+        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+  } else {
+    /* SPH gradients */
+
+    pj->primitives.gradients.rho[0] -=
+        wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+    pj->primitives.gradients.rho[1] -=
+        wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+    pj->primitives.gradients.rho[2] -=
+        wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+    pj->primitives.gradients.v[0][0] -=
+        wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+    pj->primitives.gradients.v[0][1] -=
+        wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+    pj->primitives.gradients.v[0][2] -=
+        wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+
+    pj->primitives.gradients.v[1][0] -=
+        wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+    pj->primitives.gradients.v[1][1] -=
+        wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+    pj->primitives.gradients.v[1][2] -=
+        wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+    pj->primitives.gradients.v[2][0] -=
+        wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+    pj->primitives.gradients.v[2][1] -=
+        wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+    pj->primitives.gradients.v[2][2] -=
+        wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+    pj->primitives.gradients.P[0] -=
+        wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+    pj->primitives.gradients.P[1] -=
+        wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+    pj->primitives.gradients.P[2] -=
+        wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+    pj->gravity.grad_a[0][0] -=
+        wj_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+    pj->gravity.grad_a[0][1] -=
+        wj_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+    pj->gravity.grad_a[0][2] -=
+        wj_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+
+    pj->gravity.grad_a[1][0] -=
+        wj_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+    pj->gravity.grad_a[1][1] -=
+        wj_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+    pj->gravity.grad_a[1][2] -=
+        wj_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+
+    pj->gravity.grad_a[2][0] -=
+        wj_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+    pj->gravity.grad_a[2][1] -=
+        wj_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+    pj->gravity.grad_a[2][2] -=
+        wj_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+  }
 
   hydro_slope_limit_cell_collect(pj, pi, r);
 }
@@ -250,56 +441,145 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
   xi = r * hi_inv;
   kernel_deval(xi, &wi, &wi_dx);
 
-  /* Compute gradients for pi */
-  /* there is a sign difference w.r.t. eqn. (6) because of the inverse
-   * definition of dx */
-  pi->primitives.gradients.rho[0] +=
-      (Wi[0] - Wj[0]) * wi *
-      (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-  pi->primitives.gradients.rho[1] +=
-      (Wi[0] - Wj[0]) * wi *
-      (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-  pi->primitives.gradients.rho[2] +=
-      (Wi[0] - Wj[0]) * wi *
-      (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
-  pi->primitives.gradients.v[0][0] +=
-      (Wi[1] - Wj[1]) * wi *
-      (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-  pi->primitives.gradients.v[0][1] +=
-      (Wi[1] - Wj[1]) * wi *
-      (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-  pi->primitives.gradients.v[0][2] +=
-      (Wi[1] - Wj[1]) * wi *
-      (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-  pi->primitives.gradients.v[1][0] +=
-      (Wi[2] - Wj[2]) * wi *
-      (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-  pi->primitives.gradients.v[1][1] +=
-      (Wi[2] - Wj[2]) * wi *
-      (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-  pi->primitives.gradients.v[1][2] +=
-      (Wi[2] - Wj[2]) * wi *
-      (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-  pi->primitives.gradients.v[2][0] +=
-      (Wi[3] - Wj[3]) * wi *
-      (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-  pi->primitives.gradients.v[2][1] +=
-      (Wi[3] - Wj[3]) * wi *
-      (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-  pi->primitives.gradients.v[2][2] +=
-      (Wi[3] - Wj[3]) * wi *
-      (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
-  pi->primitives.gradients.P[0] +=
-      (Wi[4] - Wj[4]) * wi *
-      (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-  pi->primitives.gradients.P[1] +=
-      (Wi[4] - Wj[4]) * wi *
-      (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-  pi->primitives.gradients.P[2] +=
-      (Wi[4] - Wj[4]) * wi *
-      (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+  if (pi->density.wcorr > const_gizmo_min_wcorr) {
+    /* Compute gradients for pi */
+    /* there is a sign difference w.r.t. eqn. (6) because of the inverse
+     * definition of dx */
+    pi->primitives.gradients.rho[0] +=
+        (Wi[0] - Wj[0]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.rho[1] +=
+        (Wi[0] - Wj[0]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.rho[2] +=
+        (Wi[0] - Wj[0]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->primitives.gradients.v[0][0] +=
+        (Wi[1] - Wj[1]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.v[0][1] +=
+        (Wi[1] - Wj[1]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.v[0][2] +=
+        (Wi[1] - Wj[1]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+    pi->primitives.gradients.v[1][0] +=
+        (Wi[2] - Wj[2]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.v[1][1] +=
+        (Wi[2] - Wj[2]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.v[1][2] +=
+        (Wi[2] - Wj[2]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+    pi->primitives.gradients.v[2][0] +=
+        (Wi[3] - Wj[3]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.v[2][1] +=
+        (Wi[3] - Wj[3]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.v[2][2] +=
+        (Wi[3] - Wj[3]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->primitives.gradients.P[0] +=
+        (Wi[4] - Wj[4]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.P[1] +=
+        (Wi[4] - Wj[4]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.P[2] +=
+        (Wi[4] - Wj[4]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->gravity.grad_a[0][0] +=
+        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->gravity.grad_a[0][1] +=
+        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->gravity.grad_a[0][2] +=
+        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->gravity.grad_a[1][0] +=
+        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->gravity.grad_a[1][1] +=
+        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->gravity.grad_a[1][2] +=
+        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->gravity.grad_a[2][0] +=
+        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->gravity.grad_a[2][1] +=
+        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->gravity.grad_a[2][2] +=
+        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+  } else {
+    /* Gradient matrix is not well-behaved, switch to SPH gradients */
+
+    pi->primitives.gradients.rho[0] -=
+        wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+    pi->primitives.gradients.rho[1] -=
+        wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+    pi->primitives.gradients.rho[2] -=
+        wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+    pi->primitives.gradients.v[0][0] -=
+        wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+    pi->primitives.gradients.v[0][1] -=
+        wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+    pi->primitives.gradients.v[0][2] -=
+        wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+    pi->primitives.gradients.v[1][0] -=
+        wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+    pi->primitives.gradients.v[1][1] -=
+        wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+    pi->primitives.gradients.v[1][2] -=
+        wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+
+    pi->primitives.gradients.v[2][0] -=
+        wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+    pi->primitives.gradients.v[2][1] -=
+        wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+    pi->primitives.gradients.v[2][2] -=
+        wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+    pi->primitives.gradients.P[0] -=
+        wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+    pi->primitives.gradients.P[1] -=
+        wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+    pi->primitives.gradients.P[2] -=
+        wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+    pi->gravity.grad_a[0][0] -=
+        wi_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+    pi->gravity.grad_a[0][1] -=
+        wi_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+    pi->gravity.grad_a[0][2] -=
+        wi_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+
+    pi->gravity.grad_a[1][0] -=
+        wi_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+    pi->gravity.grad_a[1][1] -=
+        wi_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+    pi->gravity.grad_a[1][2] -=
+        wi_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+
+    pi->gravity.grad_a[2][0] -=
+        wi_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+    pi->gravity.grad_a[2][1] -=
+        wi_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+    pi->gravity.grad_a[2][2] -=
+        wi_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+  }
 
   hydro_slope_limit_cell_collect(pi, pj, r);
 }
@@ -319,23 +599,73 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
   ih = 1.0f / h;
   const float ihdim = pow_dimension(ih);
 
-  p->primitives.gradients.rho[0] *= ihdim;
-  p->primitives.gradients.rho[1] *= ihdim;
-  p->primitives.gradients.rho[2] *= ihdim;
-
-  p->primitives.gradients.v[0][0] *= ihdim;
-  p->primitives.gradients.v[0][1] *= ihdim;
-  p->primitives.gradients.v[0][2] *= ihdim;
-  p->primitives.gradients.v[1][0] *= ihdim;
-  p->primitives.gradients.v[1][1] *= ihdim;
-  p->primitives.gradients.v[1][2] *= ihdim;
-  p->primitives.gradients.v[2][0] *= ihdim;
-  p->primitives.gradients.v[2][1] *= ihdim;
-  p->primitives.gradients.v[2][2] *= ihdim;
-
-  p->primitives.gradients.P[0] *= ihdim;
-  p->primitives.gradients.P[1] *= ihdim;
-  p->primitives.gradients.P[2] *= ihdim;
+  if (p->density.wcorr > const_gizmo_min_wcorr) {
+    p->primitives.gradients.rho[0] *= ihdim;
+    p->primitives.gradients.rho[1] *= ihdim;
+    p->primitives.gradients.rho[2] *= ihdim;
+
+    p->primitives.gradients.v[0][0] *= ihdim;
+    p->primitives.gradients.v[0][1] *= ihdim;
+    p->primitives.gradients.v[0][2] *= ihdim;
+    p->primitives.gradients.v[1][0] *= ihdim;
+    p->primitives.gradients.v[1][1] *= ihdim;
+    p->primitives.gradients.v[1][2] *= ihdim;
+    p->primitives.gradients.v[2][0] *= ihdim;
+    p->primitives.gradients.v[2][1] *= ihdim;
+    p->primitives.gradients.v[2][2] *= ihdim;
+
+    p->primitives.gradients.P[0] *= ihdim;
+    p->primitives.gradients.P[1] *= ihdim;
+    p->primitives.gradients.P[2] *= ihdim;
+
+    p->gravity.grad_a[0][0] *= ihdim;
+    p->gravity.grad_a[0][1] *= ihdim;
+    p->gravity.grad_a[0][2] *= ihdim;
+
+    p->gravity.grad_a[1][0] *= ihdim;
+    p->gravity.grad_a[1][1] *= ihdim;
+    p->gravity.grad_a[1][2] *= ihdim;
+
+    p->gravity.grad_a[2][0] *= ihdim;
+    p->gravity.grad_a[2][1] *= ihdim;
+    p->gravity.grad_a[2][2] *= ihdim;
+  } else {
+    const float ihdimp1 = pow_dimension_plus_one(ih);
+
+    float volume = p->geometry.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->primitives.gradients.v[0][0] *= ihdimp1 * volume;
+    p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
+    p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
+
+    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->primitives.gradients.v[2][0] *= ihdimp1 * volume;
+    p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
+    p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
+
+    p->primitives.gradients.P[0] *= ihdimp1 * volume;
+    p->primitives.gradients.P[1] *= ihdimp1 * volume;
+    p->primitives.gradients.P[2] *= ihdimp1 * volume;
+
+    p->gravity.grad_a[0][0] *= ihdimp1 * volume;
+    p->gravity.grad_a[0][1] *= ihdimp1 * volume;
+    p->gravity.grad_a[0][2] *= ihdimp1 * volume;
+
+    p->gravity.grad_a[1][0] *= ihdimp1 * volume;
+    p->gravity.grad_a[1][1] *= ihdimp1 * volume;
+    p->gravity.grad_a[1][2] *= ihdimp1 * volume;
+
+    p->gravity.grad_a[2][0] *= ihdimp1 * volume;
+    p->gravity.grad_a[2][1] *= ihdimp1 * volume;
+    p->gravity.grad_a[2][2] *= ihdimp1 * volume;
+  }
 
   hydro_slope_limit_cell(p);
 }
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index d707e0ee1b5707086393ea206ea9f0f60f9c1853..8798dc859a790a83ab7a3b6f1709b1302f574581 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/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);
 
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 236106a1fb04cc2e5b84f997a2389d583ce17cff..3d58be2f47c4e1904aaac5f69d1862f1d453e488 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/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);
 }
 
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index d552a3f7e86031311098293845f1aa11270c417f..6c96004847ae23b46ec3f5182f742e0e84f1118d 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/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];
 
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_face.h b/src/hydro/Gizmo/hydro_slope_limiters_face.h
index 7ae5dd2eb073d9aae8ab6f2efffdf8df15b4bb4a..ba96063d661a93a4efc4069ff7e7269a4ac58c3b 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters_face.h
+++ b/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) {
diff --git a/src/hydro/Gizmo/hydro_unphysical.h b/src/hydro/Gizmo/hydro_unphysical.h
new file mode 100644
index 0000000000000000000000000000000000000000..517e3e0918ad340580e270477c0a166590546850
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_unphysical.h
@@ -0,0 +1,55 @@
+/*******************************************************************************
+ * 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
diff --git a/src/hydro/Gizmo/hydro_velocities.h b/src/hydro/Gizmo/hydro_velocities.h
new file mode 100644
index 0000000000000000000000000000000000000000..08ba1f972b2f7a7b8a01ac4750c50a36f69784d0
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_velocities.h
@@ -0,0 +1,162 @@
+/*******************************************************************************
+ * 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;
+
+#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
+  } else {
+    /* Vacuum particles have no fluid velocity. */
+    p->v[0] = 0.;
+    p->v[1] = 0.;
+    p->v[2] = 0.;
+  }
+
+#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
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index fd05b768b232c65dafedcf0940f74d9fb9bee9e3..8f216a550ae061d01a594ff23d57575e754f85dc 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -279,7 +279,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   /* Reset the time derivatives. */
   p->u_dt = 0.0f;
   p->force.h_dt = 0.0f;
-  p->force.v_sig = p->force.soundspeed;
+  p->force.v_sig = 0.0f;
 }
 
 /**
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 90190b223a976a3283a8e589d2fcf099735473e8..4c4868cd3703e5ec5466d4878749a61284b19344 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -308,7 +308,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   p->force.h_dt = 0.0f;
 
   /* Reset maximal signal velocity */
-  p->force.v_sig = p->force.soundspeed;
+  p->force.v_sig = 0.0f;
 }
 
 /**
diff --git a/src/potential/sine_wave/potential.h b/src/potential/sine_wave/potential.h
index e2e2b8ffcc170c28a5facc8373a81746811a9991..1a4ee8aae8238c5db4c99eacb9e96bd967bcc7c4 100644
--- a/src/potential/sine_wave/potential.h
+++ b/src/potential/sine_wave/potential.h
@@ -43,6 +43,9 @@ struct external_potential {
   /*! Amplitude of the sine wave. */
   double amplitude;
 
+  /*! Growth time of the potential. */
+  double growth_time;
+
   /*! Time-step limiting factor. */
   double timestep_limit;
 };
@@ -76,7 +79,13 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
     double time, const struct external_potential* restrict potential,
     const struct phys_const* restrict phys_const, struct gpart* restrict g) {
 
-  g->a_grav[0] = potential->amplitude * sin(2. * M_PI * g->x[0]) /
+  float Acorr = 1.;
+
+  if (time < potential->growth_time) {
+    Acorr = time / potential->growth_time;
+  }
+
+  g->a_grav[0] = potential->amplitude * Acorr * sin(2. * M_PI * g->x[0]) /
                  phys_const->const_newton_G;
 }
 
@@ -114,6 +123,8 @@ static INLINE void potential_init_backend(
 
   potential->amplitude =
       parser_get_param_double(parameter_file, "SineWavePotential:amplitude");
+  potential->growth_time = parser_get_opt_param_double(
+      parameter_file, "SineWavePotential:growth_time", 0.);
   potential->timestep_limit = parser_get_param_double(
       parameter_file, "SineWavePotential:timestep_limit");
 }
diff --git a/src/riemann.h b/src/riemann.h
index 685d40708e598249151f6cbe13be016edea79553..ab6d162514326778e8d6478e07c9bae2947a7c2a 100644
--- a/src/riemann.h
+++ b/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
diff --git a/src/runner.c b/src/runner.c
index 80cb3a35c2af3448e4d3305b6bb881bd8bba3cce..54039609621945f7c529ef945c05e2ac2fe3f17c 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1918,6 +1918,8 @@ void *runner_main(void *data) {
             runner_do_recv_part(r, ci, 1, 1);
           } else if (t->subtype == task_subtype_rho) {
             runner_do_recv_part(r, ci, 1, 1);
+          } else if (t->subtype == task_subtype_gradient) {
+            runner_do_recv_part(r, ci, 1, 1);
           } else if (t->subtype == task_subtype_gpart) {
             runner_do_recv_gpart(r, ci, 1);
           } else if (t->subtype == task_subtype_spart) {
diff --git a/src/scheduler.c b/src/scheduler.c
index 7e42c3a214cff3fe30e7c885b06c48d25eac0e8d..b07c403e4ecd960b22b51f24372ca0a3420a453f 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1379,7 +1379,8 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
                           MPI_BYTE, t->ci->nodeID, t->flags, MPI_COMM_WORLD,
                           &t->req);
         } else if (t->subtype == task_subtype_xv ||
-                   t->subtype == task_subtype_rho) {
+                   t->subtype == task_subtype_rho ||
+                   t->subtype == task_subtype_gradient) {
           err = MPI_Irecv(t->ci->parts, t->ci->count, part_mpi_type,
                           t->ci->nodeID, t->flags, MPI_COMM_WORLD, &t->req);
           // message( "receiving %i parts with tag=%i from %i to %i." ,
@@ -1414,7 +1415,8 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
                           MPI_BYTE, t->cj->nodeID, t->flags, MPI_COMM_WORLD,
                           &t->req);
         } else if (t->subtype == task_subtype_xv ||
-                   t->subtype == task_subtype_rho) {
+                   t->subtype == task_subtype_rho ||
+                   t->subtype == task_subtype_gradient) {
 #ifdef SWIFT_DEBUG_CHECKS
           for (int k = 0; k < t->ci->count; k++)
             if (t->ci->parts[k].ti_drift != s->space->e->ti_current)