diff --git a/configure.ac b/configure.ac
index cbe906fa0171bdb5e139b39d6e23b8cae77244e9..95e390add27d7f45a4c6fd9a45654d73f66ebdc1 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1239,7 +1239,7 @@ esac
 # Hydro scheme.
 AC_ARG_WITH([hydro],
    [AS_HELP_STRING([--with-hydro=<scheme>],
-      [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, debug default: gadget2@:>@]
+      [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, anarchy-pu debug default: gadget2@:>@]
    )],
    [with_hydro="$withval"],
    [with_hydro="gadget2"]
@@ -1284,6 +1284,9 @@ case "$with_hydro" in
    planetary)
       AC_DEFINE([PLANETARY_SPH], [1], [Planetary SPH])
    ;;
+   anarchy-pu)
+      AC_DEFINE([ANARCHY_PU_SPH], [1], [ANARCHY (PU) SPH])
+   ;;
 
 
    *)
diff --git a/doc/RTD/source/SubgridModels/EAGLE/index.rst b/doc/RTD/source/SubgridModels/EAGLE/index.rst
index 639d98cd1a994f6f30dfc2430c90294d7486fce0..4c03a7013c21898194134a2188b638beb326366c 100644
--- a/doc/RTD/source/SubgridModels/EAGLE/index.rst
+++ b/doc/RTD/source/SubgridModels/EAGLE/index.rst
@@ -236,7 +236,7 @@ implicit problem. A valid section of the YAML file looks like:
      H_reion_z:            11.5      # Redhift of Hydrogen re-ionization
      He_reion_z_centre:     3.5      # Centre of the Gaussian used for Helium re-ionization
      He_reion_z_sigma:      0.5      # Width of the Gaussian used for Helium re-ionization
-     He_reion_eV_p_H:       2.0      # Energy injected in eV per Hydrogen atom for Helium re-ionization.
+     He_reion_ev_p_H:       2.0      # Energy injected in eV per Hydrogen atom for Helium re-ionization.
 
 And the optional parameters are:
 
diff --git a/examples/SantaBarbara/make_plots.sh b/examples/SantaBarbara/make_plots.sh
new file mode 100644
index 0000000000000000000000000000000000000000..f0ed6b3ae13acf03cb6abf8c00203462800da681
--- /dev/null
+++ b/examples/SantaBarbara/make_plots.sh
@@ -0,0 +1,4 @@
+python3 makeImage.py santabarbara_0153.hdf5 0 twilight white
+python3 plotSolution.py 153 halo
+python3 plotTempEvolution.py
+python3 rhoTHaloComparison.py
diff --git a/examples/SantaBarbara/plotSmoothingLength.py b/examples/SantaBarbara/plotSmoothingLength.py
new file mode 100644
index 0000000000000000000000000000000000000000..634cc2758d911eedbf4d2889e48011144da3d3ee
--- /dev/null
+++ b/examples/SantaBarbara/plotSmoothingLength.py
@@ -0,0 +1,162 @@
+"""
+Plots the smoothing length (compared to the softening).
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import h5py
+
+from collections import namedtuple
+
+SnapshotData = namedtuple(
+    "SnapshotData",
+    [
+        "smoothing_lengths",
+        "particle_ids",
+        "softening",
+        "internal_length",
+        "snapshot_length",
+    ],
+)
+
+HaloCatalogueData = namedtuple(
+    "HaloCatalogueData", ["largest_halo", "particle_ids_in_largest_halo"]
+)
+
+
+def load_data(filename: str) -> SnapshotData:
+    """
+    Loads the data that we need, i.e. the smoothing lengths and the
+    softening length, from the snapshot.
+    """
+
+    with h5py.File(filename, "r") as handle:
+        data = SnapshotData(
+            smoothing_lengths=handle["PartType0/SmoothingLength"][...],
+            particle_ids=handle["PartType0/ParticleIDs"][...],
+            softening=handle["GravityScheme"].attrs[
+                "Comoving softening length [internal units]"
+            ][0],
+            internal_length=handle["InternalCodeUnits"].attrs[
+                "Unit length in cgs (U_L)"
+            ][0],
+            snapshot_length=handle["Units"].attrs["Unit length in cgs (U_L)"][0],
+        )
+
+    return data
+
+
+def load_halo_data(halo_filename: str) -> HaloCatalogueData:
+    """
+    Loads the halo data and finds the particle IDs that belong to
+    the largest halo. The halo filename should be given without
+    any extension as we need a couple of files to complete this.
+    """
+
+    catalogue_filename = f"{halo_filename}.properties"
+    groups_filename = f"{halo_filename}.catalog_groups"
+    particles_filename = f"{halo_filename}.catalog_particles"
+
+    with h5py.File(catalogue_filename, "r") as handle:
+        largest_halo = np.where(
+            handle["Mass_200crit"][...] == handle["Mass_200crit"][...].max()
+        )[0][0]
+
+    with h5py.File(groups_filename, "r") as handle:
+        offset_begin = handle["Offset"][largest_halo]
+        offset_end = handle["Offset"][largest_halo + 1]
+
+    with h5py.File(particles_filename, "r") as handle:
+        particle_ids = handle["Particle_IDs"][offset_begin:offset_end]
+
+    return HaloCatalogueData(
+        largest_halo=largest_halo, particle_ids_in_largest_halo=particle_ids
+    )
+
+
+def make_plot(
+    snapshot_filename: str,
+    halo_filename: str,
+    output_filename="smoothing_length_variation.png",
+) -> None:
+    """
+    Makes the plot and saves it in output_filename.
+
+    The halo filename should be provided without extension.
+    """
+
+    data = load_data(filename)
+    halo_data = load_halo_data(halo_filename)
+
+    smoothing_lengths_in_halo = data.smoothing_lengths[
+        np.in1d(data.particle_ids, halo_data.particle_ids_in_largest_halo)
+    ]
+
+    softening = data.softening * (data.snapshot_length / data.internal_length)
+
+    fig, ax = plt.subplots(1)
+
+    ax.semilogy()
+
+    ax.hist(data.smoothing_lengths, bins="auto", label="All particles")
+    ax.hist(
+        smoothing_lengths_in_halo,
+        bins="auto",
+        label=f"Particles in largest halo (ID={halo_data.largest_halo})",
+    )
+    ax.axvline(x=softening, label="Softening", ls="--", color="grey")
+
+    ax.legend()
+
+    ax.set_xlabel("Smoothing length")
+    ax.set_ylabel("Number of particles")
+
+    ax.set_xlim(0, ax.get_xlim()[1])
+
+    fig.tight_layout()
+
+    fig.savefig(output_filename, dpi=300)
+
+    return
+
+
+if __name__ == "__main__":
+    import argparse as ap
+
+    PARSER = ap.ArgumentParser(
+        description="""
+            Makes a plot of the smoothing lengths in the box, compared
+            to the gravitational softening. Also splits out the particles
+            that are contained in the largest halo according to the
+            velociraptor outputs.
+            """
+    )
+
+    PARSER.add_argument(
+        "-s",
+        "--snapshot",
+        help="""
+            Filename and path for the snapshot (without the .hdf5),
+            Default: ./santabarbara_0153
+            """,
+        required=False,
+        default="./santabarbara_0153",
+    )
+
+    PARSER.add_argument(
+        "-v",
+        "--velociraptor",
+        help="""
+            The filename and path of the velociraptor files, excluding
+            the descriptors (i.e. without .catalog_particles).
+            Default: ./halo/santabarbara
+            """,
+        required=False,
+        default="./halo/santabarbara",
+    )
+
+    ARGS = vars(PARSER.parse_args())
+
+    filename = f"{ARGS['snapshot']}.hdf5"
+
+    make_plot(filename, ARGS["velociraptor"])
diff --git a/examples/SantaBarbara/plotSolution.py b/examples/SantaBarbara/plotSolution.py
index a23aa2089a0f82a9dad989134d1ebf11a97af2fe..131fae8f44dc3d0092e26e87a80a7862094dd4ba 100644
--- a/examples/SantaBarbara/plotSolution.py
+++ b/examples/SantaBarbara/plotSolution.py
@@ -106,13 +106,18 @@ def get_halo_data(catalogue_filename: str) -> HaloData:
     that is given by VELOCIraptor.
     """
 
+
     with h5py.File(catalogue_filename, "r") as file:
-        x = file["Xc"][0]
-        y = file["Yc"][0]
-        z = file["Zc"][0]
-        Mvir = file["Mass_200crit"][0]
-        Rvir = file["R_200crit"][0]
-        c = file["cNFW"][0]
+        largest_halo = np.where(
+            file["Mass_200crit"][...] == file["Mass_200crit"][...].max()
+        )
+
+        x = float(np.take(file["Xc"], largest_halo))
+        y = float(np.take(file["Yc"], largest_halo))
+        z = float(np.take(file["Zc"], largest_halo))
+        Mvir = float(np.take(file["Mass_200crit"], largest_halo))
+        Rvir = float(np.take(file["R_200crit"], largest_halo))
+        c = float(np.take(file["cNFW"], largest_halo))
 
     return HaloData(c=c, Rvir=Rvir, Mvir=Mvir, center=np.array([x, y, z]))
 
diff --git a/examples/SantaBarbara/plotTempEvolution.py b/examples/SantaBarbara/plotTempEvolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..90de6cb712744359dbdfbf07cc4ed81546ea38bf
--- /dev/null
+++ b/examples/SantaBarbara/plotTempEvolution.py
@@ -0,0 +1,183 @@
+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@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/>.
+#
+################################################################################
+
+# Computes the temperature evolution of the gas in a cosmological box
+
+# Physical constants needed for internal energy to temperature conversion
+k_in_J_K = 1.38064852e-23
+mH_in_kg = 1.6737236e-27
+
+# Number of snapshots generated
+n_snapshots = 153
+snapname = "santabarbara"
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+import os.path
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 9,
+'legend.fontsize': 9,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (3.15,3.15),
+'figure.subplot.left'    : 0.14,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.12,
+'figure.subplot.top'     : 0.99,
+'figure.subplot.wspace'  : 0.15,
+'figure.subplot.hspace'  : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 2.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the simulation data
+sim = h5py.File("%s_0000.hdf5" % snapname, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"][0]
+kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
+H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
+H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
+T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
+T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
+git = sim["Code"].attrs["Git Revision"]
+
+# Cosmological parameters
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
+gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
+
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
+
+unit_length_in_si = 0.01 * unit_length_in_cgs
+unit_mass_in_si = 0.001 * unit_mass_in_cgs
+unit_time_in_si = unit_time_in_cgs
+
+# Primoridal ean molecular weight as a function of temperature
+def mu(T, H_frac=H_mass_fraction, T_trans=H_transition_temp):
+    if T > T_trans:
+        return 4. / (8. - 5. * (1. - H_frac))
+    else:
+        return 4. / (1. + 3. * H_frac)
+    
+# Temperature of some primoridal gas with a given internal energy
+def T(u, H_frac=H_mass_fraction, T_trans=H_transition_temp):
+    T_over_mu = (gas_gamma - 1.) * u * mH_in_kg / k_in_J_K
+    ret = np.ones(np.size(u)) * T_trans
+
+    # Enough energy to be ionized?
+    mask_ionized = (T_over_mu > (T_trans+1) / mu(T_trans+1, H_frac, T_trans))
+    if np.sum(mask_ionized)  > 0:
+        ret[mask_ionized] = T_over_mu[mask_ionized] * mu(T_trans*10, H_frac, T_trans)
+
+    # Neutral gas?
+    mask_neutral = (T_over_mu < (T_trans-1) / mu((T_trans-1), H_frac, T_trans))
+    if np.sum(mask_neutral)  > 0:
+        ret[mask_neutral] = T_over_mu[mask_neutral] * mu(0, H_frac, T_trans)
+        
+    return ret
+
+
+z = np.zeros(n_snapshots)
+a = np.zeros(n_snapshots)
+T_mean = np.zeros(n_snapshots)
+T_std = np.zeros(n_snapshots)
+T_log_mean = np.zeros(n_snapshots)
+T_log_std = np.zeros(n_snapshots)
+T_median = np.zeros(n_snapshots)
+T_min = np.zeros(n_snapshots)
+T_max = np.zeros(n_snapshots)
+
+# Loop over all the snapshots
+for i in range(n_snapshots):
+    sim = h5py.File("%s_%04d.hdf5"% (snapname, i), "r")
+
+    z[i] = sim["/Cosmology"].attrs["Redshift"][0]
+    a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
+
+    u = sim["/PartType0/InternalEnergy"][:]
+
+    # Compute the temperature
+    u *= (unit_length_in_si**2 / unit_time_in_si**2)
+    u /= a[i]**(3 * (gas_gamma - 1.))
+    Temp = T(u)
+
+    # Gather statistics
+    T_median[i] = np.median(Temp)
+    T_mean[i] = Temp.mean()
+    T_std[i] = Temp.std()
+    T_log_mean[i] = np.log10(Temp).mean()
+    T_log_std[i] = np.log10(Temp).std()
+    T_min[i] = Temp.min()
+    T_max[i] = Temp.max()
+
+# CMB evolution
+a_evol = np.logspace(-3, 0, 60)
+T_cmb = (1. / a_evol)**2 * 2.72
+
+# Plot the interesting quantities
+figure()
+subplot(111, xscale="log", yscale="log")
+
+fill_between(a, T_mean-T_std, T_mean+T_std, color='C0', alpha=0.1)
+plot(a, T_max, ls='-.', color='C0', lw=1., label="${\\rm max}~T$")
+plot(a, T_min, ls=':', color='C0', lw=1., label="${\\rm min}~T$")
+plot(a, T_mean, color='C0', label="${\\rm mean}~T$", lw=1.5)
+fill_between(a, 10**(T_log_mean-T_log_std), 10**(T_log_mean+T_log_std), color='C1', alpha=0.1)
+plot(a, 10**T_log_mean, color='C1', label="${\\rm mean}~{\\rm log} T$", lw=1.5)
+plot(a, T_median, color='C2', label="${\\rm median}~T$", lw=1.5)
+
+legend(loc="upper left", frameon=False, handlelength=1.5)
+
+# Expected lines
+plot([1e-10, 1e10], [H_transition_temp, H_transition_temp], 'k--', lw=0.5, alpha=0.7)
+text(2.5e-2, H_transition_temp*1.07, "$T_{\\rm HII\\rightarrow HI}$", va="bottom", alpha=0.7, fontsize=8)
+plot([1e-10, 1e10], [T_minimal, T_minimal], 'k--', lw=0.5, alpha=0.7)
+text(1e-2, T_minimal*0.8, "$T_{\\rm min}$", va="top", alpha=0.7, fontsize=8)
+plot(a_evol, T_cmb, 'k--', lw=0.5, alpha=0.7)
+text(a_evol[20], T_cmb[20]*0.55, "$(1+z)^2\\times T_{\\rm CMB,0}$", rotation=-34, alpha=0.7, fontsize=8, va="top", bbox=dict(facecolor='w', edgecolor='none', pad=1.0, alpha=0.9))
+
+
+redshift_ticks = np.array([0., 1., 2., 5., 10., 20., 50., 100.])
+redshift_labels = ["$0$", "$1$", "$2$", "$5$", "$10$", "$20$", "$50$", "$100$"]
+a_ticks = 1. / (redshift_ticks + 1.)
+
+xticks(a_ticks, redshift_labels)
+minorticks_off()
+
+xlabel("${\\rm Redshift}~z$", labelpad=0)
+ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
+xlim(9e-3, 1.1)
+ylim(20, 2.5e7)
+
+savefig("Temperature_evolution.png", dpi=200)
+
diff --git a/examples/SantaBarbara/rhoTHaloComparison.py b/examples/SantaBarbara/rhoTHaloComparison.py
new file mode 100644
index 0000000000000000000000000000000000000000..edf4e47527bd61e3f9b017b5a53510036dbcacac
--- /dev/null
+++ b/examples/SantaBarbara/rhoTHaloComparison.py
@@ -0,0 +1,221 @@
+"""
+This script finds the temperatures inside all of the halos and
+compares it against the virial temperature. This uses velociraptor
+and the SWIFT snapshot.
+
+Folkert Nobels (2018) nobels@strw.leidenuniv.nl
+Josh Borrow (2019) joshua.borrow@durham.ac.uk
+"""
+
+import numpy as np
+import h5py
+import matplotlib.pyplot as plt
+from matplotlib.colors import LogNorm
+
+mH = 1.6733e-24  # g
+kB = 1.38e-16  # erg/K
+
+
+def virial_temp(mu, M, h=0.703, a=1.0):
+    """
+    Calculates the virial temperature according to
+
+    https://arxiv.org/pdf/1105.5701.pdf
+
+    Equation 1.
+    """
+    return  4e4 * (mu / 1.2) * (M * h / 1e8) ** (2 / 3) / (10 * a)
+
+
+def calculate_group_sizes_array(offsets: np.array, total_size: int) -> np.array:
+    """
+    Calculates the group sizes array from the offsets and total size, i.e. it
+    calculates the diff between all of the offsets.
+    """
+
+    # Does not include the LAST one
+    group_sizes = [x - y for x, y in zip(offsets[1:], offsets[:-1])]
+    group_sizes += [total_size - offsets[-1]]
+    group_sizes = np.array(group_sizes, dtype=type(offsets[0]))
+
+    return group_sizes
+
+
+def create_group_array(group_sizes: np.array) -> np.array:
+    """
+    Creates an array that looks like:
+    [GroupID0, GroupID0, ..., GroupIDN, GroupIDN]
+    i.e. for each group create the correct number of group ids.
+    This is used to be sorted alongside the particle IDs to track
+    the placement of group IDs.
+    """
+
+    slices = []
+    running_total_of_particles = type(offsets[0])(0)
+
+    for group in group_sizes:
+        slices.append([running_total_of_particles, group + running_total_of_particles])
+        running_total_of_particles += group
+
+    groups = np.empty(group_sizes.sum(), dtype=int)
+
+    for group_id, group in enumerate(slices):
+        groups[group[0] : group[1]] = group_id
+
+    return groups
+
+
+if __name__ == "__main__":
+    import argparse as ap
+
+    PARSER = ap.ArgumentParser(
+        description="""
+        Makes many plots comparing the virial temperature and the
+        temperature of halos. Requires the velociraptor files and
+        the SWIFT snapshot.
+        """
+    )
+
+    PARSER.add_argument(
+        "-s",
+        "--snapshot",
+        help="""
+            Filename and path for the snapshot (without the .hdf5),
+            Default: ./santabarbara_0153
+            """,
+        required=False,
+        default="./santabarbara_0153",
+    )
+
+    PARSER.add_argument(
+        "-v",
+        "--velociraptor",
+        help="""
+            The filename and path of the velociraptor files, excluding
+            the descriptors (i.e. without .catalog_particles).
+            Default: ./halo/santabarbara
+            """,
+        required=False,
+        default="./halo/santabarbara",
+    )
+
+    ARGS = vars(PARSER.parse_args())
+
+    # Grab some metadata before we begin.
+    with h5py.File("%s.hdf5" % ARGS["snapshot"], "r") as handle:
+        # Cosmology
+        a = handle["Cosmology"].attrs["Scale-factor"][0]
+        h = handle["Cosmology"].attrs["h"][0]
+
+        # Gas
+        hydro = handle["HydroScheme"].attrs
+        X = hydro["Hydrogen mass fraction"][0]
+        gamma = hydro["Adiabatic index"][0]
+        mu = 1 / (X + (1 - X) / 4)
+
+    # First we must construct a group array so we know which particles belong
+    # to which group.
+    with h5py.File("%s.catalog_groups" % ARGS["velociraptor"], "r") as handle:
+        offsets = handle["Offset"][...]
+
+    # Then, extract the particles that belong to the halos. For that, we need
+    # the particle IDs:
+    with h5py.File("%s.catalog_particles" % ARGS["velociraptor"], "r") as handle:
+        ids_in_halos = handle["Particle_IDs"][...]
+
+    number_of_groups = len(offsets)
+    group_sizes = calculate_group_sizes_array(offsets, ids_in_halos.size)
+    group_array = create_group_array(group_sizes)
+
+    # We can now load the particle data from the snapshot.
+    with h5py.File("%s.hdf5" % ARGS["snapshot"], "r") as handle:
+        gas_particles = handle["PartType0"]
+
+        particle_ids = gas_particles["ParticleIDs"][...]
+
+        # Requires numpy 1.15 or greater.
+        _, particles_in_halos_mask, group_array_mask = np.intersect1d(
+            particle_ids,
+            ids_in_halos,
+            assume_unique=True,
+            return_indices=True,
+        )
+
+        # We also need to re-index the group array to cut out DM particles
+        group_array = group_array[group_array_mask]
+        
+        # Kill the spare
+        del particle_ids
+
+        # Now we can only read the properties that we require from the snapshot!
+        temperatures = np.take(gas_particles["InternalEnergy"], particles_in_halos_mask)
+        # This 1e10 should probably be explained somewhere...
+        temperatures *= 1e10 * (gamma - 1) * mu * mH / kB
+
+        densities = np.take(gas_particles["Density"], particles_in_halos_mask)
+
+    # Just a quick check to make sure nothing's gone wrong.
+    assert len(group_array) == len(temperatures)
+
+    # Now we can loop through all the particles and find out the mean temperature and
+    # density in each halo.
+
+    particles_in_group = np.zeros(number_of_groups, dtype=int)
+    temp_in_group = np.zeros(number_of_groups, dtype=float)
+    dens_in_group = np.zeros(number_of_groups, dtype=float)
+
+    for group, T, rho in zip(group_array, temperatures, densities):
+        particles_in_group[group] += 1
+        temp_in_group[group] += T
+        dens_in_group[group] += rho
+
+    # First get a mask to ensure no runtime warnings
+    mask = particles_in_group != 0
+    
+    # Normalize
+    temp_in_group[mask] /= particles_in_group[mask]
+    dens_in_group[mask] /= particles_in_group[mask]
+
+    # Now we can load the data according to the halo finder to compare with.
+    with h5py.File("%s.properties" % ARGS["velociraptor"], "r") as handle:
+        halo_masses = handle["Mass_200crit"][...]
+
+    halo_temperatures = virial_temp(mu, halo_masses * 1e10, h=h, a=a)
+
+    # Finally, the plotting!
+
+    fig, ax = plt.subplots()
+    ax.loglog()
+
+    mask = np.logical_and.reduce([
+         halo_temperatures != 0.0,
+         temp_in_group != 0.0,
+    ])
+
+    temp_in_group = temp_in_group[mask]
+    halo_temperatures = halo_temperatures[mask]
+
+    mean_range = [temp_in_group.min(), temp_in_group.max()]
+    halo_range = [halo_temperatures.min(), halo_temperatures.max()]
+
+    bottom = min([halo_range[0], mean_range[0]])
+    top = max([halo_range[1], mean_range[1]])
+
+    plt.plot(
+        [bottom, top],
+        [bottom, top],
+        lw=2, linestyle="--", color="grey", label="1:1"
+    )
+    
+    ax.scatter(halo_temperatures, temp_in_group, s=2, edgecolor="none", label="Halos")
+
+    ax.set_ylabel("Mean Group Temperature [K]")
+    ax.set_xlabel("Halo Virial Temperature [K]")
+
+    ax.set_ylim(mean_range)
+    ax.set_xlim(halo_range)
+
+    ax.legend(frameon=False)
+
+    fig.tight_layout()
+    fig.savefig("temperature_comparison.png", dpi=300)
diff --git a/examples/SantaBarbara/velociraptor_cfg.cfg b/examples/SantaBarbara/velociraptor_cfg.cfg
index 4b9b441b71cc65a85a9731018d38ffc2f003c0ff..b904d2a419ad6010e12073209bd77e8f59eef7c4 100644
--- a/examples/SantaBarbara/velociraptor_cfg.cfg
+++ b/examples/SantaBarbara/velociraptor_cfg.cfg
@@ -11,7 +11,19 @@ Cosmological_input=1
 #Type of snapshot to read. Ignored when using within SWIFT.
 HDF_name_convention=6 # ILLUSTRIS 0, GADGETX 1, EAGLE 2, GIZMO 3, SIMBA 4, MUFASA 5, SWIFTEAGLE 6
 
-Particle_search_type=2 #search all particles, see allvars for other types
+#whether star particles are present in the input
+Input_includes_star_particle=0
+#bhs present
+Input_includes_bh_particle=0
+#no wind present
+Input_includes_wind_particle=0
+#no tracers present
+Input_includes_tracer_particle=0
+#no low res/extra dm particle types present
+Input_includes_extradm_particle=0
+
+
+Particle_search_type=1 #search all particles, see allvars for other types
 Baryon_searchflag=0 #if 1 search for baryons separately using phase-space search when identifying substructures, 2 allows special treatment in field FOF linking and phase-space substructure search, 0 treat the same as dark matter particles
 Search_for_substructure=0 #if 0, end search once field objects are found
 FoF_Field_search_type=5 #5 3DFOF search for field halos, 4 for 6DFOF clean up of field halos, 3 for 6DFOF with velocity scale distinct for each halo
diff --git a/examples/SedovBlast_1D/plotSolution.py b/examples/SedovBlast_1D/plotSolution.py
index 2738b7c8f301a7351d962ac0f29faccd0a770fc9..c6d4a989da7493f7b500946610eea8832696bf6f 100644
--- a/examples/SedovBlast_1D/plotSolution.py
+++ b/examples/SedovBlast_1D/plotSolution.py
@@ -83,6 +83,12 @@ S = sim["/PartType0/Entropy"][:]
 P = sim["/PartType0/Pressure"][:]
 rho = sim["/PartType0/Density"][:]
 
+try:
+    alpha = sim["/PartType0/Viscosity"][:]
+    plot_alpha = True 
+except:
+    plot_alpha = False
+
 
 # Now, work our the solution....
 
@@ -246,14 +252,23 @@ ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
 xlim(0, 1.3 * r_shock)
 ylim(-2, 22)
 
-# Entropy profile ---------------------------------
+# Entropy/alpha profile ---------------------------------
 subplot(235)
-plot(r, S, '.', color='r', ms=2.)
-plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+
+if plot_alpha:
+    plot(r, alpha, '.', color='r', ms=2.0)
+    plot([r_shock, r_shock], [-1, 2], "--", color="k", alpha=0.8, lw=1.2)
+    ylabel(r"${\rm{Viscosity}}~\alpha$", labelpad=0)
+    # Show location of shock
+    ylim(0, 2)
+else:
+    plot(r, S, '.', color='r', ms=2.0)
+    plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+    ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+    ylim(-5, 50)
+
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=0)
 xlim(0, 1.3 * r_shock)
-ylim(-5, 50)
 
 # Information -------------------------------------
 subplot(236, frameon=False)
diff --git a/examples/SodShock_1D/plotSolution.py b/examples/SodShock_1D/plotSolution.py
index 0149d66a0c28c777a4265da10d86ed160086995d..12ae8a9cf26fb715281ed00bf565c5c8d9a234fb 100644
--- a/examples/SodShock_1D/plotSolution.py
+++ b/examples/SodShock_1D/plotSolution.py
@@ -87,6 +87,11 @@ try:
     plot_alpha = True 
 except:
     plot_alpha = False
+try:
+    alpha_diff = sim["PartType0/Diffusion"][:]
+    plot_alpha_diff = True
+except:
+    plot_alpha_diff = False
 
 N = 1000  # Number of points
 x_min = -1.
@@ -239,12 +244,20 @@ ylim(-0.1, 0.95)
 
 # Density profile --------------------------------
 subplot(232)
-plot(x, rho, '.', color='r', ms=4.0)
-plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+if plot_alpha_diff:
+    plot(x, alpha_diff, '.', color='r', ms=4.0)
+    ylabel(r"${\rm{Diffusion}}~\alpha$", labelpad=0)
+    # Show location of contact discontinuity
+    plot([x_34, x_34], [-100, 100], color="k", alpha=0.5, ls="dashed", lw=1.2)
+    ylim(0, 1)
+else:
+    plot(x, rho, '.', color='r', ms=4.0)
+    plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+    ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+    ylim(0.05, 1.1)
+
 xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
 xlim(-0.5, 0.5)
-ylim(0.05, 1.1)
 
 # Pressure profile --------------------------------
 subplot(233)
@@ -257,7 +270,7 @@ ylim(0.01, 1.1)
 
 # Internal energy profile -------------------------
 subplot(234)
-plot(x, u, '.', color='r', ms=4.0)
+scatter(x, u, marker='.', c=alpha_diff, s=4.0)
 plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
 xlabel("${\\rm{Position}}~x$", labelpad=0)
 ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
diff --git a/examples/SodShock_1D/sodShock.yml b/examples/SodShock_1D/sodShock.yml
index 69554b4db733166fc5dbb6d198966fd8f9b8d49c..b936f50f6c2c7d9078c49fbad868aa5334498957 100644
--- a/examples/SodShock_1D/sodShock.yml
+++ b/examples/SodShock_1D/sodShock.yml
@@ -27,6 +27,10 @@ Statistics:
 SPH:
   resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
   CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  viscosity_alpha_min:   0.01
+  viscosity_alpha:       0.01
+  viscosity_alpha_max:   2.0
+  viscosity_length:      0.02
 
 # Parameters related to the initial conditions
 InitialConditions:
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 75f7d77ac8659d8c0d8a75f622691a4ca21ae772..759fbe4aa60e42c434bb6b75ebf0e7b0fd531ff0 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -35,9 +35,14 @@ SPH:
   H_mass_fraction:       0.755    # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy. Default value is derived from the physical constants.
   H_ionization_temperature: 1e4   # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
   viscosity_alpha:       0.8      # (Optional) Override for the initial value of the artificial viscosity. In schemes that have a fixed AV, this remains as alpha throughout the run.
-  viscosity_alpha_max:   2.0      # (Optional) Maximal value for the artificial viscosity in schemes that allow alpha to vary
-  viscosity_alpha_min:   0.1      # (Optional) Minimal value for the artificial viscosity in schemes that allow alpha to vary
-  viscosity_length:      0.1      # (Optional) Decay length for the artificial viscosity in schemes that allow alpha to vary
+  viscosity_alpha_max:   2.0      # (Optional) Maximal value for the artificial viscosity in schemes that allow alpha to vary.
+  viscosity_alpha_min:   0.1      # (Optional) Minimal value for the artificial viscosity in schemes that allow alpha to vary.
+  viscosity_length:      0.1      # (Optional) Decay length for the artificial viscosity in schemes that allow alpha to vary.
+  diffusion_alpha:       0.0      # (Optional) Override the initial value for the thermal diffusion coefficient in schemes with thermal diffusion.
+  diffusion_beta:        0.01     # (Optional) Override the decay/rise rate tuning parameter for the thermal diffusion.
+  diffusion_alpha_max:   1.0      # (Optional) Override the maximal thermal diffusion coefficient that is allowed for a given particle.
+  diffusion_alpha_min:   0.0      # (Optional) Override the minimal thermal diffusion coefficient that is allowed for a given particle.
+
 
 # Parameters for the self-gravity scheme
 Gravity:
diff --git a/src/debug.c b/src/debug.c
index 809d7048c45888eb41bd277bdb4971d469a98dc3..2c1d81b676520c4823a27f3a7a8cc617585ca5f2 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -60,6 +60,8 @@
 #include "./hydro/Shadowswift/hydro_debug.h"
 #elif defined(PLANETARY_SPH)
 #include "./hydro/Planetary/hydro_debug.h"
+#elif defined(ANARCHY_PU_SPH)
+#include "./hydro/AnarchyPU/hydro_debug.h"
 #else
 #error "Invalid choice of SPH variant"
 #endif
diff --git a/src/engine.c b/src/engine.c
index 5d51d3b7a32d15c33a01c01d1810fabbd78bc6c2..fd76580722e4481b716d6fd6ec010cdfcaa07a5a 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2643,7 +2643,8 @@ void engine_skip_force_and_kick(struct engine *e) {
         t->subtype == task_subtype_limiter || t->subtype == task_subtype_grav ||
         t->type == task_type_end_force ||
         t->type == task_type_grav_long_range || t->type == task_type_grav_mm ||
-        t->type == task_type_grav_down || t->type == task_type_cooling)
+        t->type == task_type_grav_down || t->type == task_type_cooling ||
+        t->type == task_type_extra_ghost || t->subtype == task_subtype_gradient)
       t->skip = 1;
   }
 
diff --git a/src/hydro.h b/src/hydro.h
index 15c45c1dcfa1217d842904dfc1303acea607e3ab..3bf7d2228b528796a8717d6a5ab17fda6f569d25 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -72,6 +72,11 @@
 #include "./hydro/Planetary/hydro.h"
 #include "./hydro/Planetary/hydro_iact.h"
 #define SPH_IMPLEMENTATION "Minimal version of SPH with multiple materials"
+#elif defined(ANARCHY_PU_SPH)
+#include "./hydro/AnarchyPU/hydro.h"
+#include "./hydro/AnarchyPU/hydro_iact.h"
+#define SPH_IMPLEMENTATION \
+  "ANARCHY (Pressure-Energy) SPH (Dalla Vecchia+ in prep)"
 #else
 #error "Invalid choice of SPH variant"
 #endif
diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h
new file mode 100644
index 0000000000000000000000000000000000000000..27501944c8afd6cd316717158057086a8f6ffd4e
--- /dev/null
+++ b/src/hydro/AnarchyPU/hydro.h
@@ -0,0 +1,930 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
+ *                    Josh Borrow (joshua.borrow@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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_ANARCHY_PU_HYDRO_H
+#define SWIFT_ANARCHY_PU_HYDRO_H
+
+/**
+ * @file PressureEnergy/hydro.h
+ * @brief P-U conservative implementation of SPH (Non-neighbour loop
+ * equations)
+ *
+ * The thermal variable is the internal energy (u). A simple constant
+ * viscosity term with a Balsara switch is implemented.
+ *
+ * No thermal conduction term is implemented.
+ *
+ * This implementation corresponds to the one presented in the SWIFT
+ * documentation and in Hopkins, "A general class of Lagrangian smoothed
+ * particle hydrodynamics methods and implications for fluid mixing problems",
+ * MNRAS, 2013.
+ */
+
+#include "adiabatic_index.h"
+#include "approx_math.h"
+#include "cosmology.h"
+#include "dimension.h"
+#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "hydro_space.h"
+#include "kernel_hydro.h"
+#include "minmax.h"
+
+#include <float.h>
+
+/**
+ * @brief Returns the comoving internal energy of a particle at the last
+ * time the particle was kicked.
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy(const struct part *restrict p,
+                                   const struct xpart *restrict xp) {
+
+  return xp->u_full;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle at the last
+ * time the particle was kicked.
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy(const struct part *restrict p,
+                                   const struct xpart *restrict xp,
+                                   const struct cosmology *cosmo) {
+
+  return xp->u_full * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the comoving internal energy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_internal_energy(const struct part *restrict p) {
+
+  return p->u;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_physical_internal_energy(const struct part *restrict p,
+                                           const struct cosmology *cosmo) {
+
+  return p->u * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
+    const struct part *restrict p) {
+
+  return p->pressure_bar;
+}
+
+/**
+ * @brief Returns the physical pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties and
+ * convert it to physical coordinates.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
+    const struct part *restrict p, const struct cosmology *cosmo) {
+
+  return cosmo->a_factor_pressure * p->pressure_bar;
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle at the last
+ * time the particle was kicked.
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
+    const struct part *restrict p, const struct xpart *restrict xp) {
+
+  return gas_entropy_from_internal_energy(p->rho, xp->u_full);
+}
+
+/**
+ * @brief Returns the physical entropy of a particle at the last
+ * time the particle was kicked.
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
+    const struct part *restrict p, const struct xpart *restrict xp,
+    const struct cosmology *cosmo) {
+
+  /* Note: no cosmological conversion required here with our choice of
+   * coordinates. */
+  return gas_entropy_from_internal_energy(p->rho, xp->u_full);
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_entropy(const struct part *restrict p) {
+
+  return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the physical entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_physical_entropy(const struct part *restrict p,
+                                   const struct cosmology *cosmo) {
+
+  /* Note: no cosmological conversion required here with our choice of
+   * coordinates. */
+  return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the comoving sound speed of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part *restrict p) {
+
+  /* Compute the sound speed -- see theory section for justification */
+  /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
+  const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho);
+
+  return square_rooted;
+}
+
+/**
+ * @brief Returns the physical sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_soundspeed(const struct part *restrict p,
+                              const struct cosmology *cosmo) {
+
+  return cosmo->a_factor_sound_speed * p->force.soundspeed;
+}
+
+/**
+ * @brief Returns the comoving density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
+    const struct part *restrict p) {
+
+  return p->rho;
+}
+
+/**
+ * @brief Returns the comoving density of a particle.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
+    const struct part *restrict p, const struct cosmology *cosmo) {
+
+  return cosmo->a3_inv * p->rho;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+    const struct part *restrict p) {
+
+  return p->mass;
+}
+
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+    struct part *restrict p, float m) {
+
+  p->mass = m;
+}
+
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
+ * @param dt_kick_grav The time (for gravity accelerations) since the last kick.
+ * @param v (return) The velocities at the current time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
+    const struct part *restrict p, const struct xpart *xp, float dt_kick_hydro,
+    float dt_kick_grav, float v[3]) {
+
+  v[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro +
+         xp->a_grav[0] * dt_kick_grav;
+  v[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro +
+         xp->a_grav[1] * dt_kick_grav;
+  v[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro +
+         xp->a_grav[2] * dt_kick_grav;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy_dt(const struct part *restrict p) {
+
+  return p->u_dt;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest
+ * @param cosmo Cosmology data structure
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy_dt(const struct part *restrict p,
+                                      const struct cosmology *cosmo) {
+
+  return p->u_dt * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Sets the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest.
+ * @param du_dt The new time derivative of the internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_comoving_internal_energy_dt(struct part *restrict p, float du_dt) {
+
+  p->u_dt = du_dt;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest.
+ * @param cosmo Cosmology data structure
+ * @param du_dt The new time derivative of the internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_physical_internal_energy_dt(struct part *restrict p,
+                                      const struct cosmology *cosmo,
+                                      float du_dt) {
+
+  p->u_dt = du_dt / cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Computes the hydro time-step of a given particle
+ *
+ * This function returns the time-step of a particle given its hydro-dynamical
+ * state. A typical time-step calculation would be the use of the CFL condition.
+ *
+ * @param p Pointer to the particle data
+ * @param xp Pointer to the extended particle data
+ * @param hydro_properties The SPH parameters
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
+    const struct part *restrict p, const struct xpart *restrict xp,
+    const struct hydro_props *restrict hydro_properties,
+    const struct cosmology *restrict cosmo) {
+
+  const float CFL_condition = hydro_properties->CFL_condition;
+
+  /* CFL condition */
+  const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
+                       (cosmo->a_factor_sound_speed * p->viscosity.v_sig);
+
+  const float dt_u_change =
+      (p->u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->u_dt) : FLT_MAX;
+
+  return fminf(dt_cfl, dt_u_change);
+}
+
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+    struct part *p, float dt) {}
+
+/**
+ * @brief Operations performed when a particle gets removed from the
+ * simulation volume.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void hydro_remove_part(
+    const struct part *p, const struct xpart *xp) {}
+
+/**
+ * @brief Prepares a particle for the density calculation.
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the various density loop over neighbours. Typically, all fields of the
+ * density sub-structure of a particle get zeroed in here.
+ *
+ * @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
+ */
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+    struct part *restrict p, const struct hydro_space *hs) {
+
+  p->density.wcount = 0.f;
+  p->density.wcount_dh = 0.f;
+  p->rho = 0.f;
+  p->density.rho_dh = 0.f;
+  p->pressure_bar = 0.f;
+  p->density.pressure_bar_dh = 0.f;
+
+  p->density.rot_v[0] = 0.f;
+  p->density.rot_v[1] = 0.f;
+  p->density.rot_v[2] = 0.f;
+
+  p->viscosity.div_v = 0.f;
+  p->diffusion.laplace_u = 0.f;
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * Multiplies the density and number of neighbours by the appropiate constants
+ * and add the self-contribution term.
+ * Additional quantities such as velocity gradients will also get the final
+ * terms added to them here.
+ *
+ * Also adds/multiplies the cosmological terms if need be.
+ *
+ * @param p The particle to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+    struct part *restrict p, const struct cosmology *cosmo) {
+
+  /* Some smoothing length multiples. */
+  const float h = p->h;
+  const float h_inv = 1.0f / h;                       /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
+  const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+  /* Final operation on the density (add self-contribution). */
+  p->rho += p->mass * kernel_root;
+  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
+  p->pressure_bar += p->mass * p->u * kernel_root;
+  p->density.pressure_bar_dh -= hydro_dimension * p->mass * p->u * kernel_root;
+  p->density.wcount += kernel_root;
+  p->density.wcount_dh -= hydro_dimension * kernel_root;
+
+  /* Finish the calculation by inserting the missing h-factors */
+  p->rho *= h_inv_dim;
+  p->density.rho_dh *= h_inv_dim_plus_one;
+  p->pressure_bar *= (h_inv_dim * hydro_gamma_minus_one);
+  p->density.pressure_bar_dh *= (h_inv_dim_plus_one * hydro_gamma_minus_one);
+  p->density.wcount *= h_inv_dim;
+  p->density.wcount_dh *= h_inv_dim_plus_one;
+
+  const float rho_inv = 1.f / p->rho;
+  const float a_inv2 = cosmo->a2_inv;
+
+  /* Finish calculation of the velocity curl components */
+  p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+  p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+  p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+
+  /* Finish calculation of the velocity divergence */
+  p->viscosity.div_v *=
+      h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension;
+}
+
+/**
+ * @brief Prepare a particle for the gradient calculation.
+ *
+ * This function is called after the density loop and before the gradient loop.
+ *
+ * We use it to set the physical timestep for the particle and to copy the
+ * actual velocities, which we need to boost our interfaces during the flux
+ * calculation. We also initialize the variables used for the time step
+ * calculation.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
+    struct part *restrict p, struct xpart *restrict xp,
+    const struct cosmology *cosmo) {
+
+  const float fac_B = cosmo->a_factor_Balsara_eps;
+
+  /* Compute the norm of the curl */
+  const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
+                             p->density.rot_v[1] * p->density.rot_v[1] +
+                             p->density.rot_v[2] * p->density.rot_v[2]);
+
+  /* Compute the norm of div v */
+  const float abs_div_v = fabsf(p->viscosity.div_v);
+
+  /* Compute the sound speed -- see theory section for justification */
+  const float soundspeed = hydro_get_comoving_soundspeed(p);
+
+  /* Compute the Balsara switch */
+  const float balsara =
+      abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_B / p->h);
+
+  /* Compute the "grad h" term */
+  const float common_factor = p->h / (hydro_dimension * p->density.wcount);
+  const float grad_h_term = (p->density.pressure_bar_dh * common_factor *
+                             hydro_one_over_gamma_minus_one) /
+                            (1.f + common_factor * p->density.wcount_dh);
+
+  /* Update variables. */
+  p->force.f = grad_h_term;
+  p->force.soundspeed = soundspeed;
+  p->force.balsara = balsara;
+}
+
+/**
+ * @brief Resets the variables that are required for a gradient calculation.
+ *
+ * This function is called after hydro_prepare_gradient.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_gradient(
+    struct part *restrict p) {
+  p->viscosity.v_sig = 2.f * p->force.soundspeed;
+}
+
+/**
+ * @brief Finishes the gradient calculation.
+ *
+ * Just a wrapper around hydro_gradients_finalize, which can be an empty method,
+ * in which case no gradients are used.
+ *
+ * This method also initializes the force loop variables.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_gradient(
+    struct part *p) {
+
+  /* Some smoothing length multiples. */
+  const float h = p->h;
+  const float h_inv = 1.0f / h;                       /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
+  const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+  /* Include the extra factors in the del^2 u */
+
+  p->diffusion.laplace_u *= 2 * h_inv_dim_plus_one;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * In the desperate case where a particle has no neighbours (likely because
+ * of the h_max ceiling), set the particle fields to something sensible to avoid
+ * NaNs in the next calculations.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
+    struct part *restrict p, struct xpart *restrict xp,
+    const struct cosmology *cosmo) {
+
+  /* Some smoothing length multiples. */
+  const float h = p->h;
+  const float h_inv = 1.0f / h;                 /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+  /* Re-set problematic values */
+  p->rho = p->mass * kernel_root * h_inv_dim;
+  p->viscosity.v_sig = 0.f;
+  p->pressure_bar =
+      p->mass * p->u * hydro_gamma_minus_one * kernel_root * h_inv_dim;
+  p->density.wcount = kernel_root * h_inv_dim;
+  p->density.rho_dh = 0.f;
+  p->density.wcount_dh = 0.f;
+  p->density.pressure_bar_dh = 0.f;
+
+  p->density.rot_v[0] = 0.f;
+  p->density.rot_v[1] = 0.f;
+  p->density.rot_v[2] = 0.f;
+
+  /* Probably not shocking, so this is safe to do */
+  p->viscosity.div_v = 0.f;
+  p->diffusion.laplace_u = 0.f;
+}
+
+/**
+ * @brief Prepare a particle for the force calculation.
+ *
+ * This function is called in the ghost task to convert some quantities coming
+ * from the density loop over neighbours into quantities ready to be used in the
+ * force loop over neighbours. Quantities are typically read from the density
+ * sub-structure and written to the force sub-structure.
+ * Examples of calculations done here include the calculation of viscosity term
+ * constants, thermal conduction terms, hydro conversions, etc.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties.
+ * @param dt_alpha The time-step used to evolve non-cosmological quantities such
+ *                 as the artificial viscosity.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+    struct part *restrict p, struct xpart *restrict xp,
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
+
+  /* Here we need to update the artificial viscosity */
+
+  /* Timescale for decay */
+  const float tau =
+      p->h / (2.f * p->viscosity.v_sig * hydro_props->viscosity.length);
+  /* Construct time differential of div.v implicitly */
+  const float div_v_dt =
+      dt_alpha == 0.f
+          ? 0.f
+          : (p->viscosity.div_v - p->viscosity.div_v_previous_step) / dt_alpha;
+  /* Construct the source term for the AV; if shock detected this is _positive_
+   * as div_v_dt should be _negative_ before the shock hits */
+  const float S = p->h * p->h * max(0.f, -1.f * div_v_dt);
+  const float v_sig_square = p->viscosity.v_sig * p->viscosity.v_sig;
+  /* Calculate the current appropriate value of the AV based on the above */
+  const float alpha_loc =
+      hydro_props->viscosity.alpha_max * S / (v_sig_square + S);
+
+  if (alpha_loc > p->viscosity.alpha) {
+    /* Reset the value of alpha to the appropriate value */
+    p->viscosity.alpha = alpha_loc;
+  } else {
+    /* Integrate the alpha forward in time to decay back to alpha = 0 */
+    const float alpha_dt = (alpha_loc - p->viscosity.alpha) / tau;
+
+    /* Finally, we can update the actual value of the alpha */
+    p->viscosity.alpha += alpha_dt * dt_alpha;
+  }
+
+  if (p->viscosity.alpha < hydro_props->viscosity.alpha_min) {
+    p->viscosity.alpha = hydro_props->viscosity.alpha_min;
+  }
+
+  /* Set our old div_v to the one for the next loop */
+  p->viscosity.div_v_previous_step = p->viscosity.div_v;
+
+  /* Now for the diffusive alpha */
+
+  const float sqrt_u = sqrtf(p->u);
+  /* Calculate initial value of alpha dt before bounding */
+  /* alpha_diff_dt is cosmology-less */
+  float alpha_diff_dt =
+      hydro_props->diffusion.beta * p->h * p->diffusion.laplace_u / sqrt_u;
+
+  float new_diffusion_alpha = p->diffusion.alpha + alpha_diff_dt * dt_alpha;
+
+  if (new_diffusion_alpha > hydro_props->diffusion.alpha_max) {
+    new_diffusion_alpha = hydro_props->diffusion.alpha_max;
+  } else if (new_diffusion_alpha < hydro_props->diffusion.alpha_min) {
+    new_diffusion_alpha = hydro_props->diffusion.alpha_min;
+  }
+
+  p->diffusion.alpha = new_diffusion_alpha;
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * Resets all hydro acceleration and time derivative fields in preparation
+ * for the sums taking  place in the various force tasks.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+    struct part *restrict p) {
+
+  /* Reset the acceleration. */
+  p->a_hydro[0] = 0.0f;
+  p->a_hydro[1] = 0.0f;
+  p->a_hydro[2] = 0.0f;
+
+  /* Reset the time derivatives. */
+  p->u_dt = 0.0f;
+  p->force.h_dt = 0.0f;
+}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+    struct part *restrict p, const struct xpart *restrict xp) {
+
+  /* Re-set the predicted velocities */
+  p->v[0] = xp->v_full[0];
+  p->v[1] = xp->v_full[1];
+  p->v[2] = xp->v_full[2];
+
+  /* Re-set the entropy */
+  p->u = xp->u_full;
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * Additional hydrodynamic quantites are drifted forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * Note the different time-step sizes used for the different quantities as they
+ * include cosmological factors.
+ *
+ * @param p The particle.
+ * @param xp The extended data of the particle.
+ * @param dt_drift The drift time-step for positions.
+ * @param dt_therm The drift time-step for thermal quantities.
+ */
+__attribute__((always_inline)) INLINE static void hydro_predict_extra(
+    struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
+    float dt_therm) {
+
+  const float h_inv = 1.f / p->h;
+
+  /* Predict smoothing length */
+  const float w1 = p->force.h_dt * h_inv * dt_drift;
+  if (fabsf(w1) < 0.2f)
+    p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+  else
+    p->h *= expf(w1);
+
+  /* Predict density and weighted pressure */
+  const float w2 = -hydro_dimension * w1;
+  if (fabsf(w2) < 0.2f) {
+    const float expf_approx =
+        approx_expf(w2); /* 4th order expansion of exp(w) */
+    p->rho *= expf_approx;
+    p->pressure_bar *= expf_approx;
+  } else {
+    const float expf_exact = expf(w2);
+    p->rho *= expf_exact;
+    p->pressure_bar *= expf_exact;
+  }
+
+  /* Predict the internal energy */
+  p->u += p->u_dt * dt_therm;
+
+  /* Compute the new sound speed */
+  const float soundspeed = hydro_get_comoving_soundspeed(p);
+
+  p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Finishes the force calculation.
+ *
+ * Multiplies the force and accelerations by the appropiate constants
+ * and add the self-contribution term. In most cases, there is little
+ * to do here.
+ *
+ * Cosmological terms are also added/multiplied here.
+ *
+ * @param p The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+    struct part *restrict p, const struct cosmology *cosmo) {
+
+  p->force.h_dt *= p->h * hydro_dimension_inv;
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * Additional hydrodynamic quantites are kicked forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * @param p The particle to act upon.
+ * @param xp The particle extended data to act upon.
+ * @param dt_therm The time-step for this kick (for thermodynamic quantities).
+ * @param dt_grav The time-step for this kick (for gravity quantities).
+ * @param dt_hydro The time-step for this kick (for hydro quantities).
+ * @param dt_kick_corr The time-step for this kick (for gravity corrections).
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+    struct part *restrict p, struct xpart *restrict xp, float dt_therm,
+    float dt_grav, float dt_hydro, float dt_kick_corr,
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
+
+  /* Do not decrease the energy by more than a factor of 2*/
+  if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * xp->u_full) {
+    p->u_dt = -0.5f * xp->u_full / dt_therm;
+  }
+  xp->u_full += p->u_dt * dt_therm;
+
+  /* Apply the minimal energy limit */
+  const float min_energy =
+      hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+  if (xp->u_full < min_energy) {
+    xp->u_full = min_energy;
+    p->u_dt = 0.f;
+  }
+
+  /* Compute the sound speed */
+  const float soundspeed = hydro_get_comoving_soundspeed(p);
+
+  p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Converts hydro quantity of a particle at the start of a run
+ *
+ * This function is called once at the end of the engine_init_particle()
+ * routine (at the start of a calculation) after the densities of
+ * particles have been computed.
+ * This can be used to convert internal energy into entropy for instance.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle to act upon
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme.
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+    struct part *restrict p, struct xpart *restrict xp,
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
+
+  /* Convert the physcial internal energy to the comoving one. */
+  /* u' = a^(3(g-1)) u */
+  const float factor = 1.f / cosmo->a_factor_internal_energy;
+  p->u *= factor;
+  xp->u_full = p->u;
+
+  /* Apply the minimal energy limit */
+  const float min_energy =
+      hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+  if (xp->u_full < min_energy) {
+    xp->u_full = min_energy;
+    p->u = min_energy;
+    p->u_dt = 0.f;
+  }
+
+  /* Note that unlike Minimal the pressure and sound speed cannot be calculated
+   * here because they are smoothed properties in this scheme. */
+
+  /* Set the initial value of the artificial viscosity based on the non-variable
+     schemes for safety */
+
+  p->viscosity.alpha = hydro_props->viscosity.alpha;
+  /* Initialise this here to keep all the AV variables together */
+  p->viscosity.div_v_previous_step = 0.f;
+
+  /* Set the initial values for the thermal diffusion */
+  p->diffusion.alpha = hydro_props->diffusion.alpha;
+}
+
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions or assignments between the particle
+ * and extended particle fields.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+    struct part *restrict p, struct xpart *restrict xp) {
+
+  p->time_bin = 0;
+  p->wakeup = time_bin_not_awake;
+  xp->v_full[0] = p->v[0];
+  xp->v_full[1] = p->v[1];
+  xp->v_full[2] = p->v[2];
+  xp->a_grav[0] = 0.f;
+  xp->a_grav[1] = 0.f;
+  xp->a_grav[2] = 0.f;
+  xp->u_full = p->u;
+
+  hydro_reset_acceleration(p);
+  hydro_init_part(p, NULL);
+}
+
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+  p->u = u_init;
+}
+
+#endif /* SWIFT_ANARCHY_PU_HYDRO_H */
diff --git a/src/hydro/AnarchyPU/hydro_debug.h b/src/hydro/AnarchyPU/hydro_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..79ab5b96653a3f503c1baf255f4296f0ccc4aca9
--- /dev/null
+++ b/src/hydro/AnarchyPU/hydro_debug.h
@@ -0,0 +1,44 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
+ *                    Josh Borrow (joshua.borrow@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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_ANARCHY_PU_HYDRO_DEBUG_H
+#define SWIFT_ANARCHY_PU_HYDRO_DEBUG_H
+/**
+ * @file PressureEnergy/hydro_debug.h
+ * @brief P-U conservative implementation of SPH (Debugging routines)
+ */
+
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+    const struct part* p, const struct xpart* xp) {
+  printf(
+      "x=[%.3e,%.3e,%.3e], "
+      "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
+      "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
+      "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, \n"
+      "p_dh=%.3e, p_bar=%.3e, alpha=%.3e \n"
+      "time_bin=%d\n",
+      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
+      xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+      p->u, p->u_dt, p->viscosity.v_sig, hydro_get_comoving_pressure(p), p->h,
+      p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
+      p->density.pressure_bar_dh, p->pressure_bar, p->viscosity.alpha,
+      p->time_bin);
+}
+
+#endif /* SWIFT_ANARCHY_PU_HYDRO_DEBUG_H */
diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..c214db3b018e00b7f3881fb301b55d6cf49a1f43
--- /dev/null
+++ b/src/hydro/AnarchyPU/hydro_iact.h
@@ -0,0 +1,577 @@
+/*******************************************************************************
+ * This file is part* of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
+ *                    Josh Borrow (joshua.borrow@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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_ANARCHY_PU_HYDRO_IACT_H
+#define SWIFT_ANARCHY_PU_HYDRO_IACT_H
+
+/**
+ * @file PressureEnergy/hydro_iact.h
+ * @brief P-U implementation of SPH (Neighbour loop equations)
+ *
+ * The thermal variable is the internal energy (u). A simple constant
+ * viscosity term with a Balsara switch is implemented.
+ *
+ * No thermal conduction term is implemented.
+ *
+ * See PressureEnergy/hydro.h for references.
+ */
+
+#include "adiabatic_index.h"
+#include "minmax.h"
+
+/**
+ * @brief Density interaction between two part*icles.
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+    float r2, const float* dx, float hi, float hj, struct part* pi,
+    struct part* pj, float a, float H) {
+
+  float wi, wj, wi_dx, wj_dx;
+  float dv[3], curlvr[3];
+
+  const float r = sqrtf(r2);
+
+  /* Get the masses. */
+  const float mi = pi->mass;
+  const float mj = pj->mass;
+
+  /* Compute density of pi. */
+  const float hi_inv = 1.f / hi;
+  const float ui = r * hi_inv;
+
+  kernel_deval(ui, &wi, &wi_dx);
+
+  pi->rho += mj * wi;
+  pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+
+  pi->pressure_bar += mj * wi * pj->u;
+  pi->density.pressure_bar_dh -=
+      mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
+  pi->density.wcount += wi;
+  pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+  /* Compute density of pj. */
+  const float hj_inv = 1.f / hj;
+  const float uj = r * hj_inv;
+  kernel_deval(uj, &wj, &wj_dx);
+
+  pj->rho += mi * wj;
+  pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+  pj->pressure_bar += mi * wj * pi->u;
+  pj->density.pressure_bar_dh -=
+      mi * pi->u * (hydro_dimension * wj + uj * wj_dx);
+  pj->density.wcount += wj;
+  pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+
+  /* Now we need to compute the div terms */
+  const float r_inv = 1.f / r;
+  const float faci = mj * wi_dx * r_inv;
+  const float facj = mi * wj_dx * r_inv;
+
+  /* Compute dv dot r */
+  dv[0] = pi->v[0] - pj->v[0];
+  dv[1] = pi->v[1] - pj->v[1];
+  dv[2] = pi->v[2] - pj->v[2];
+  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+  pi->viscosity.div_v -= faci * dvdr;
+  pj->viscosity.div_v -= facj * dvdr;
+
+  /* Compute dv cross r */
+  curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+  curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+  curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+  pi->density.rot_v[0] += faci * curlvr[0];
+  pi->density.rot_v[1] += faci * curlvr[1];
+  pi->density.rot_v[2] += faci * curlvr[2];
+
+  /* Negative because of the change in sign of dx & dv. */
+  pj->density.rot_v[0] += facj * curlvr[0];
+  pj->density.rot_v[1] += facj * curlvr[1];
+  pj->density.rot_v[2] += facj * curlvr[2];
+}
+
+/**
+ * @brief Density interaction between two part*icles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+    float r2, const float* dx, float hi, float hj, struct part* pi,
+    const struct part* pj, float a, float H) {
+
+  float wi, wi_dx;
+  float dv[3], curlvr[3];
+
+  /* Get the masses. */
+  const float mj = pj->mass;
+
+  /* Get r and r inverse. */
+  const float r = sqrtf(r2);
+
+  const float h_inv = 1.f / hi;
+  const float ui = r * h_inv;
+  kernel_deval(ui, &wi, &wi_dx);
+
+  pi->rho += mj * wi;
+  pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+
+  pi->pressure_bar += mj * wi * pj->u;
+
+  pi->density.pressure_bar_dh -=
+      mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
+  pi->density.wcount += wi;
+  pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+  const float r_inv = 1.f / r;
+  const float faci = mj * wi_dx * r_inv;
+
+  /* Compute dv dot r */
+  dv[0] = pi->v[0] - pj->v[0];
+  dv[1] = pi->v[1] - pj->v[1];
+  dv[2] = pi->v[2] - pj->v[2];
+  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+  pi->viscosity.div_v -= faci * dvdr;
+
+  /* Compute dv cross r */
+  curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+  curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+  curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+  pi->density.rot_v[0] += faci * curlvr[0];
+  pi->density.rot_v[1] += faci * curlvr[1];
+  pi->density.rot_v[2] += faci * curlvr[2];
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j
+ *
+ * This method wraps around hydro_gradients_collect, which can be an empty
+ * method, in which case no gradients are used.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_gradient(
+    float r2, const float* dx, float hi, float hj, struct part* restrict pi,
+    struct part* restrict pj, float a, float H) {
+
+  /* We need to construct the maximal signal velocity between our particle
+   * and all of it's neighbours */
+
+  const float dv_dx = (pi->v[0] - pj->v[0]) * dx[0] +
+                      (pi->v[1] - pj->v[1]) * dx[1] +
+                      (pi->v[2] - pj->v[2]) * dx[2];
+
+  const float dv_dx_factor = min(0, const_viscosity_beta * dv_dx);
+
+  const float new_v_sig =
+      pi->force.soundspeed + pj->force.soundspeed - dv_dx_factor;
+
+  /* Update if we need to */
+  pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig);
+  pj->viscosity.v_sig = max(pj->viscosity.v_sig, new_v_sig);
+
+  /* Calculate Del^2 u for the thermal diffusion coefficient. */
+  /* Need to get some kernel values F_ij = wi_dx */
+  float wi, wi_dx, wj, wj_dx;
+
+  const float r = sqrtf(r2);
+  const float ui = r / hi;
+  const float uj = r / hj;
+
+  kernel_deval(ui, &wi, &wi_dx);
+  kernel_deval(uj, &wj, &wj_dx);
+
+  const float delta_u_factor = (pi->u - pj->u) / r;
+  pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho;
+  pj->diffusion.laplace_u -= pi->mass * delta_u_factor * wj_dx / pi->rho;
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j:
+ * non-symmetric version
+ *
+ * This method wraps around hydro_gradients_nonsym_collect, which can be an
+ * empty method, in which case no gradients are used.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
+    float r2, const float* dx, float hi, float hj, struct part* restrict pi,
+    struct part* restrict pj, float a, float H) {
+
+  /* We need to construct the maximal signal velocity between our particle
+   * and all of it's neighbours */
+
+  const float dv_dx = (pi->v[0] - pj->v[0]) * dx[0] +
+                      (pi->v[1] - pj->v[1]) * dx[1] +
+                      (pi->v[2] - pj->v[2]) * dx[2];
+
+  const float dv_dx_factor = min(0, const_viscosity_beta * dv_dx);
+
+  const float new_v_sig =
+      pi->force.soundspeed + pj->force.soundspeed - dv_dx_factor;
+
+  /* Update if we need to */
+  pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig);
+
+  /* Calculate Del^2 u for the thermal diffusion coefficient. */
+  /* Need to get some kernel values F_ij = wi_dx */
+  float wi, wi_dx;
+
+  const float r = sqrtf(r2);
+  const float ui = r / hi;
+
+  kernel_deval(ui, &wi, &wi_dx);
+
+  const float delta_u_factor = (pi->u - pj->u) / r;
+  pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho;
+}
+
+/**
+ * @brief Force interaction between two part*icles.
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+    float r2, const float* dx, float hi, float hj, struct part* pi,
+    struct part* pj, float a, float H) {
+
+  /* Cosmological factors entering the EoMs */
+  const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+  const float a2_Hubble = a * a * H;
+
+  const float r = sqrtf(r2);
+  const float r_inv = 1.0f / r;
+
+  /* Recover some data */
+  const float mj = pj->mass;
+  const float mi = pi->mass;
+
+  const float miui = mi * pi->u;
+  const float mjuj = mj * pj->u;
+
+  const float rhoi = pi->rho;
+  const float rhoj = pj->rho;
+  /* Compute gradient terms */
+  const float f_ij = 1.f - (pi->force.f / mjuj);
+  const float f_ji = 1.f - (pj->force.f / miui);
+
+  /* Get the kernel for hi. */
+  const float hi_inv = 1.0f / hi;
+  const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+  const float xi = r * hi_inv;
+  float wi, wi_dx;
+  kernel_deval(xi, &wi, &wi_dx);
+  const float wi_dr = hid_inv * wi_dx;
+
+  /* Get the kernel for hj. */
+  const float hj_inv = 1.0f / hj;
+  const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+  const float xj = r * hj_inv;
+  float wj, wj_dx;
+  kernel_deval(xj, &wj, &wj_dx);
+  const float wj_dr = hjd_inv * wj_dx;
+
+  /* Compute dv dot r. */
+  const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+                     (pi->v[1] - pj->v[1]) * dx[1] +
+                     (pi->v[2] - pj->v[2]) * dx[2];
+
+  /* Includes the hubble flow term; not used for du/dt */
+  const float dvdr_Hubble = dvdr + a2_Hubble * r2;
+
+  /* Are the part*icles moving towards each others ? */
+  const float omega_ij = min(dvdr_Hubble, 0.f);
+  const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+  /* Compute sound speeds and signal velocity */
+  const float v_sig = 0.5 * (pi->viscosity.v_sig + pj->viscosity.v_sig);
+
+  /* Balsara term */
+  const float balsara_i = pi->force.balsara;
+  const float balsara_j = pj->force.balsara;
+
+  /* Construct the full viscosity term */
+  const float rho_ij = rhoi + rhoj;
+  const float alpha = pi->viscosity.alpha + pj->viscosity.alpha;
+  const float visc =
+      -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
+
+  /* Convolve with the kernel */
+  const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+
+  /* SPH acceleration term */
+  const float sph_acc_term =
+      pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one *
+      ((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) *
+      r_inv;
+
+  /* Assemble the acceleration */
+  const float acc = sph_acc_term + visc_acc_term;
+
+  /* Use the force Luke ! */
+  pi->a_hydro[0] -= mj * acc * dx[0];
+  pi->a_hydro[1] -= mj * acc * dx[1];
+  pi->a_hydro[2] -= mj * acc * dx[2];
+
+  pj->a_hydro[0] += mi * acc * dx[0];
+  pj->a_hydro[1] += mi * acc * dx[1];
+  pj->a_hydro[2] += mi * acc * dx[2];
+
+  /* Get the time derivative for u. */
+  const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
+                              pj->u * pi->u * (f_ij / pi->pressure_bar) *
+                              wi_dr * dvdr * r_inv;
+  const float sph_du_term_j = hydro_gamma_minus_one * hydro_gamma_minus_one *
+                              pi->u * pj->u * (f_ji / pj->pressure_bar) *
+                              wj_dr * dvdr * r_inv;
+
+  /* Viscosity term */
+  const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
+
+  /* Diffusion term */
+  const float v_diff =
+      max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
+  const float alpha_diff = 0.5 * (pi->diffusion.alpha + pj->diffusion.alpha);
+  /* wi_dx + wj_dx / 2 is F_ij */
+  const float diff_du_term =
+      alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij;
+
+  /* Assemble the energy equation term */
+  const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
+  const float du_dt_j = sph_du_term_j + visc_du_term - diff_du_term;
+
+  /* Internal energy time derivative */
+  pi->u_dt += du_dt_i * mj;
+  pj->u_dt += du_dt_j * mi;
+
+  /* Get the time derivative for h. */
+  pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+  pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
+}
+
+/**
+ * @brief Force interaction between two part*icles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+    float r2, const float* dx, float hi, float hj, struct part* pi,
+    const struct part* pj, float a, float H) {
+
+  /* Cosmological factors entering the EoMs */
+  const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+  const float a2_Hubble = a * a * H;
+
+  const float r = sqrtf(r2);
+  const float r_inv = 1.0f / r;
+
+  /* Recover some data */
+  // const float mi = pi->mass;
+  const float mj = pj->mass;
+  const float mi = pi->mass;
+
+  const float miui = mi * pi->u;
+  const float mjuj = mj * pj->u;
+
+  const float rhoi = pi->rho;
+  const float rhoj = pj->rho;
+  /* Compute gradient terms */
+  const float f_ij = 1.f - (pi->force.f / mjuj);
+  const float f_ji = 1.f - (pj->force.f / miui);
+
+  /* Get the kernel for hi. */
+  const float hi_inv = 1.0f / hi;
+  const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+  const float xi = r * hi_inv;
+  float wi, wi_dx;
+  kernel_deval(xi, &wi, &wi_dx);
+  const float wi_dr = hid_inv * wi_dx;
+
+  /* Get the kernel for hj. */
+  const float hj_inv = 1.0f / hj;
+  const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+  const float xj = r * hj_inv;
+  float wj, wj_dx;
+  kernel_deval(xj, &wj, &wj_dx);
+  const float wj_dr = hjd_inv * wj_dx;
+
+  /* Compute dv dot r. */
+  const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+                     (pi->v[1] - pj->v[1]) * dx[1] +
+                     (pi->v[2] - pj->v[2]) * dx[2];
+
+  /* Includes the hubble flow term; not used for du/dt */
+  const float dvdr_Hubble = dvdr + a2_Hubble * r2;
+
+  /* Are the part*icles moving towards each others ? */
+  const float omega_ij = min(dvdr_Hubble, 0.f);
+  const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+  /* Compute sound speeds and signal velocity */
+  const float v_sig = 0.5 * (pi->viscosity.v_sig + pj->viscosity.v_sig);
+
+  /* Balsara term */
+  const float balsara_i = pi->force.balsara;
+  const float balsara_j = pj->force.balsara;
+
+  /* Construct the full viscosity term */
+  const float rho_ij = rhoi + rhoj;
+  const float alpha = pi->viscosity.alpha + pj->viscosity.alpha;
+  const float visc =
+      -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
+
+  /* Convolve with the kernel */
+  const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+
+  /* SPH acceleration term */
+  const float sph_acc_term =
+      pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one *
+      ((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) *
+      r_inv;
+
+  /* Assemble the acceleration */
+  const float acc = sph_acc_term + visc_acc_term;
+
+  /* Use the force Luke ! */
+  pi->a_hydro[0] -= mj * acc * dx[0];
+  pi->a_hydro[1] -= mj * acc * dx[1];
+  pi->a_hydro[2] -= mj * acc * dx[2];
+
+  /* Get the time derivative for u. */
+  const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
+                              pj->u * pi->u * (f_ij / pi->pressure_bar) *
+                              wi_dr * dvdr * r_inv;
+
+  /* Viscosity term */
+  const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
+
+  /* Diffusion term */
+  const float v_diff =
+      max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
+  const float alpha_diff = 0.5 * (pi->diffusion.alpha + pj->diffusion.alpha);
+  /* wi_dx + wj_dx / 2 is F_ij */
+  const float diff_du_term =
+      alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij;
+
+  /* Assemble the energy equation term */
+  const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
+
+  /* Internal energy time derivative */
+  pi->u_dt += du_dt_i * mj;
+
+  /* Get the time derivative for h. */
+  pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+}
+
+/**
+ * @brief Timestep limiter loop
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_limiter(
+    float r2, const float* dx, float hi, float hj, struct part* restrict pi,
+    struct part* restrict pj, float a, float H) {
+
+  /* Nothing to do here if both particles are active */
+}
+
+/**
+ * @brief Timestep limiter loop (non-symmetric version)
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_limiter(
+    float r2, const float* dx, float hi, float hj, struct part* restrict pi,
+    struct part* restrict pj, float a, float H) {
+
+  /* Wake up the neighbour? */
+  if (pi->viscosity.v_sig >
+      const_limiter_max_v_sig_ratio * pj->viscosity.v_sig) {
+
+    pj->wakeup = time_bin_awake;
+  }
+}
+
+#endif /* SWIFT_MINIMAL_HYDRO_IACT_H */
diff --git a/src/hydro/AnarchyPU/hydro_io.h b/src/hydro/AnarchyPU/hydro_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..99fea6a9feb2722da3c82482d4f5e79330e5c7e6
--- /dev/null
+++ b/src/hydro/AnarchyPU/hydro_io.h
@@ -0,0 +1,224 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
+ *                    Josh Borrow (joshua.borrow@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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_ANARCHY_PU_HYDRO_IO_H
+#define SWIFT_ANARCHY_PU_HYDRO_IO_H
+/**
+ * @file PressureEnergy/hydro_io.h
+ * @brief P-U implementation of SPH (i/o routines)
+ *
+ * The thermal variable is the internal energy (u). A simple constant
+ * viscosity term with a Balsara switch is implemented.
+ *
+ * No thermal conduction term is implemented.
+ *
+ * See PressureEnergy/hydro.h for references.
+ */
+
+#include "adiabatic_index.h"
+#include "hydro.h"
+#include "io_properties.h"
+#include "kernel_hydro.h"
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void hydro_read_particles(struct part* parts,
+                                        struct io_props* list,
+                                        int* num_fields) {
+
+  *num_fields = 8;
+
+  /* List what we want to read */
+  list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+                                UNIT_CONV_LENGTH, parts, x);
+  list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+                                UNIT_CONV_SPEED, parts, v);
+  list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+                                parts, mass);
+  list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+                                UNIT_CONV_LENGTH, parts, h);
+  list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+                                UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+  list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+                                UNIT_CONV_NO_UNITS, parts, id);
+  list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
+                                UNIT_CONV_ACCELERATION, parts, a_hydro);
+  list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
+                                UNIT_CONV_DENSITY, parts, rho);
+}
+
+INLINE static void convert_S(const struct engine* e, const struct part* p,
+                             const struct xpart* xp, float* ret) {
+
+  ret[0] = hydro_get_comoving_entropy(p, xp);
+}
+
+INLINE static void convert_P(const struct engine* e, const struct part* p,
+                             const struct xpart* xp, float* ret) {
+
+  ret[0] = hydro_get_comoving_pressure(p);
+}
+
+INLINE static void convert_part_pos(const struct engine* e,
+                                    const struct part* p,
+                                    const struct xpart* xp, double* ret) {
+
+  if (e->s->periodic) {
+    ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
+    ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
+    ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
+  } else {
+    ret[0] = p->x[0];
+    ret[1] = p->x[1];
+    ret[2] = p->x[2];
+  }
+}
+
+INLINE static void convert_part_vel(const struct engine* e,
+                                    const struct part* p,
+                                    const struct xpart* xp, float* ret) {
+
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+  const struct cosmology* cosmo = e->cosmology;
+  const integertime_t ti_current = e->ti_current;
+  const double time_base = e->time_base;
+
+  const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+  const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+  /* Get time-step since the last kick */
+  float dt_kick_grav, dt_kick_hydro;
+  if (with_cosmology) {
+    dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+    dt_kick_grav -=
+        cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+    dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+    dt_kick_hydro -=
+        cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+  } else {
+    dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+    dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+  }
+
+  /* Extrapolate the velocites to the current time */
+  hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+  /* Conversion from internal units to peculiar velocities */
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
+}
+
+INLINE static void convert_part_potential(const struct engine* e,
+                                          const struct part* p,
+                                          const struct xpart* xp, float* ret) {
+  if (p->gpart != NULL)
+    ret[0] = gravity_get_comoving_potential(p->gpart);
+  else
+    ret[0] = 0.f;
+}
+
+INLINE static void convert_viscosity(const struct engine* e,
+                                     const struct part* p,
+                                     const struct xpart* xp, float* ret) {
+  ret[0] = p->viscosity.alpha;
+}
+
+INLINE static void convert_diffusion(const struct engine* e,
+                                     const struct part* p,
+                                     const struct xpart* xp, float* ret) {
+  ret[0] = p->diffusion.alpha;
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+INLINE static void hydro_write_particles(const struct part* parts,
+                                         const struct xpart* xparts,
+                                         struct io_props* list,
+                                         int* num_fields) {
+
+  *num_fields = 12;
+
+  /* List what we want to write */
+  list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+                                              UNIT_CONV_LENGTH, parts, xparts,
+                                              convert_part_pos);
+  list[1] = io_make_output_field_convert_part(
+      "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
+  list[2] =
+      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
+  list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+                                 parts, h);
+  list[4] = io_make_output_field("InternalEnergy", FLOAT, 1,
+                                 UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+  list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+                                 UNIT_CONV_NO_UNITS, parts, id);
+  list[6] =
+      io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+  list[7] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
+                                 parts, pressure_bar);
+  list[8] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+                                              UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+                                              parts, xparts, convert_S);
+  list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+                                              UNIT_CONV_POTENTIAL, parts,
+                                              xparts, convert_part_potential);
+  list[10] = io_make_output_field_convert_part("Viscosity", FLOAT, 1,
+                                               UNIT_CONV_NO_UNITS, parts,
+                                               xparts, convert_viscosity);
+  list[11] = io_make_output_field_convert_part("Diffusion", FLOAT, 1,
+                                               UNIT_CONV_NO_UNITS, parts,
+                                               xparts, convert_diffusion);
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+INLINE static void hydro_write_flavour(hid_t h_grpsph) {
+
+  /* Viscosity and thermal conduction */
+  /* Nothing in this minimal model... */
+  io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
+  io_write_attribute_s(h_grpsph, "Viscosity Model",
+                       "Minimal treatment as in Monaghan (1992)");
+
+  /* Time integration properties */
+  io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
+                       const_max_u_change);
+}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+INLINE static int writeEntropyFlag(void) { return 0; }
+
+#endif /* SWIFT_ANARCHY_PU_HYDRO_IO_H */
diff --git a/src/hydro/AnarchyPU/hydro_part.h b/src/hydro/AnarchyPU/hydro_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..e25ca4426376576157390e7da6d5bb69b22c9418
--- /dev/null
+++ b/src/hydro/AnarchyPU/hydro_part.h
@@ -0,0 +1,210 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
+ *                    Josh Borrow (joshua.borrow@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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_ANARCHY_PU_HYDRO_PART_H
+#define SWIFT_ANARCHY_PU_HYDRO_PART_H
+/**
+ * @file PressureEnergy/hydro_part.h
+ * @brief P-U implementation of SPH (Particle definition)
+ *
+ * The thermal variable is the internal energy (u). A simple constant
+ * viscosity term with a Balsara switch is implemented.
+ *
+ * No thermal conduction term is implemented.
+ *
+ * See PressureEnergy/hydro.h for references.
+ */
+
+#include "chemistry_struct.h"
+#include "cooling_struct.h"
+
+/**
+ * @brief Particle fields not needed during the SPH loops over neighbours.
+ *
+ * This structure contains the particle fields that are not used in the
+ * density or force loops. Quantities should be used in the kick, drift and
+ * potentially ghost tasks only.
+ */
+struct xpart {
+
+  /*! Offset between current position and position at last tree rebuild. */
+  float x_diff[3];
+
+  /*! Offset between the current position and position at the last sort. */
+  float x_diff_sort[3];
+
+  /*! Velocity at the last full step. */
+  float v_full[3];
+
+  /*! Gravitational acceleration at the last full step. */
+  float a_grav[3];
+
+  /*! Internal energy at the last full step. */
+  float u_full;
+
+  /*! Additional data used to record cooling information */
+  struct cooling_xpart_data cooling_data;
+
+} SWIFT_STRUCT_ALIGN;
+
+/**
+ * @brief Particle fields for the SPH particles
+ *
+ * The density and force substructures are used to contain variables only used
+ * within the density and force loops over neighbours. All more permanent
+ * variables should be declared in the main part of the part structure,
+ */
+struct part {
+
+  /*! Particle unique ID. */
+  long long id;
+
+  /*! Pointer to corresponding gravity part. */
+  struct gpart* gpart;
+
+  /*! Particle position. */
+  double x[3];
+
+  /*! Particle predicted velocity. */
+  float v[3];
+
+  /*! Particle acceleration. */
+  float a_hydro[3];
+
+  /*! Particle mass. */
+  float mass;
+
+  /*! Particle smoothing length. */
+  float h;
+
+  /*! Particle internal energy. */
+  float u;
+
+  /*! Time derivative of the internal energy. */
+  float u_dt;
+
+  /*! Particle density. */
+  float rho;
+
+  /*! Particle pressure (weighted) */
+  float pressure_bar;
+
+  /* Store viscosity information in a separate struct. */
+  struct {
+
+    /*! Particle velocity divergence */
+    float div_v;
+
+    /*! Particle velocity divergence from previous step */
+    float div_v_previous_step;
+
+    /*! Artificial viscosity parameter */
+    float alpha;
+
+    /*! Signal velocity */
+    float v_sig;
+
+  } viscosity;
+
+  /* Store thermal diffusion information in a separate struct. */
+  struct {
+
+    /*! del^2 u, a smoothed quantity */
+    float laplace_u;
+
+    /*! Thermal diffusion coefficient */
+    float alpha;
+
+  } diffusion;
+
+  /* Store density/force specific stuff. */
+  union {
+
+    /**
+     * @brief Structure for the variables only used in the density loop over
+     * neighbours.
+     *
+     * Quantities in this sub-structure should only be accessed in the density
+     * loop over neighbours and the ghost task.
+     */
+    struct {
+
+      /*! Neighbour number count. */
+      float wcount;
+
+      /*! Derivative of the neighbour number with respect to h. */
+      float wcount_dh;
+
+      /*! Derivative of density with respect to h */
+      float rho_dh;
+
+      /*! Derivative of the weighted pressure with respect to h */
+      float pressure_bar_dh;
+
+      /*! Particle velocity curl. */
+      float rot_v[3];
+
+    } density;
+
+    /**
+     * @brief Structure for the variables only used in the force loop over
+     * neighbours.
+     *
+     * Quantities in this sub-structure should only be accessed in the force
+     * loop over neighbours and the ghost, drift and kick tasks.
+     */
+    struct {
+
+      /*! "Grad h" term -- only partial in P-U */
+      float f;
+
+      /*! Particle soundspeed. */
+      float soundspeed;
+
+      /*! Time derivative of smoothing length  */
+      float h_dt;
+
+      /*! Balsara switch */
+      float balsara;
+
+    } force;
+  };
+
+  /* Chemistry information */
+  struct chemistry_part_data chemistry_data;
+
+  /*! Time-step length */
+  timebin_t time_bin;
+
+  /* Need waking up ? */
+  char wakeup;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+  /* Time of the last drift */
+  integertime_t ti_drift;
+
+  /* Time of the last kick */
+  integertime_t ti_kick;
+
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_ANARCHY_PU_HYDRO_PART_H */
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 524774435d03a6d808c4535a6c54b68ad16bcb66..2af6eb47e8de091f391dc758e7842cc5cada702a 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -334,7 +334,7 @@ hydro_set_physical_internal_energy_dt(struct part *restrict p,
                                       const struct cosmology *cosmo,
                                       float du_dt) {
 
-  p->u_dt = du_dt * cosmo->a_factor_internal_energy;
+  p->u_dt = du_dt / cosmo->a_factor_internal_energy;
 }
 
 /**
diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h
index ed7aa6b89d2b50ab2e00cedb0b3ef6779689feb1..68dfa3eb825b7dd984eb6fd52336ee97b0c4f93e 100644
--- a/src/hydro/Planetary/hydro.h
+++ b/src/hydro/Planetary/hydro.h
@@ -357,7 +357,7 @@ hydro_set_physical_internal_energy_dt(struct part *restrict p,
                                       const struct cosmology *cosmo,
                                       float du_dt) {
 
-  p->u_dt = du_dt * cosmo->a_factor_internal_energy;
+  p->u_dt = du_dt / cosmo->a_factor_internal_energy;
 }
 
 /**
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index 400a84915b700464b9b86f74400ba578b4efa446..7b6c302fb7f8001cb96d595eb35fb6ef08222bec 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -359,7 +359,7 @@ hydro_set_physical_internal_energy_dt(struct part *restrict p,
                                       const struct cosmology *cosmo,
                                       float du_dt) {
 
-  p->u_dt = du_dt * cosmo->a_factor_internal_energy;
+  p->u_dt = du_dt / cosmo->a_factor_internal_energy;
 }
 
 /**
@@ -535,7 +535,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     const struct cosmology *cosmo, const struct hydro_props *hydro_props,
     const float dt_alpha) {
 
-  const float fac_mu = cosmo->a_factor_mu;
+  const float fac_B = cosmo->a_factor_Balsara_eps;
 
   /* Compute the norm of the curl */
   const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
@@ -551,7 +551,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   /* Compute the Balsara switch */
   const float balsara =
       hydro_props->viscosity.alpha * abs_div_v /
-      (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu / p->h);
+      (abs_div_v + curl_v + 0.0001f * soundspeed * fac_B / p->h);
 
   /* Compute the "grad h" term */
   const float common_factor = p->h / (hydro_dimension * p->density.wcount);
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
index 7ef55b86c24972f8f287273441da99f26285c531..9767a165d3f9243c9dd9b72bdc9131fe338e52b4 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
@@ -498,7 +498,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     const struct cosmology *cosmo, const struct hydro_props *hydro_props,
     const float dt_alpha) {
 
-  const float fac_mu = cosmo->a_factor_mu;
+  const float fac_B = cosmo->a_factor_Balsara_eps;
 
   const float h_inv = 1.f / p->h;
 
@@ -515,7 +515,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Compute the Balsara switch */
   const float balsara =
-      abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu * h_inv);
+      abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_B * h_inv);
 
   /* Compute the "grad h" term */
   const float common_factor = p->h / (hydro_dimension * p->density.wcount);
diff --git a/src/hydro_io.h b/src/hydro_io.h
index 1a2d6319b7caf6c09b9af406cbdd323f27607791..60e5593cc0630ef4bc33ab407f6a669b7de1def1 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -43,6 +43,8 @@
 #include "./hydro/Shadowswift/hydro_io.h"
 #elif defined(PLANETARY_SPH)
 #include "./hydro/Planetary/hydro_io.h"
+#elif defined(ANARCHY_PU_SPH)
+#include "./hydro/AnarchyPU/hydro_io.h"
 #else
 #error "Invalid choice of SPH variant"
 #endif
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 167c13a30c22b01c20e355d1b8c60903ca026ad8..e83b90cdd932e2a61baaace72922025cf6e199f6 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -41,12 +41,30 @@
 #define hydro_props_default_min_temp 0.f
 #define hydro_props_default_H_ionization_temperature 1e4
 #define hydro_props_default_viscosity_alpha 0.8f
+
+#ifdef ANARCHY_PU_SPH
+/* This nasty #ifdef is only temporary until we separate the viscosity
+ * and hydro components. If it is not removed by July 2019, shout at JB. */
 #define hydro_props_default_viscosity_alpha_min \
-  0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
+  0.01f /* values taken from Schaller+ 2015 */
+#define hydro_props_default_viscosity_alpha_max \
+  2.0f /* values taken from Schaller+ 2015 */
+#define hydro_props_default_viscosity_length \
+  0.01f /* values taken from Schaller+ 2015 */
+#else
+#define hydro_props_default_viscosity_alpha_min \
+  0.1f /* values taken from (price,2004), not used in legacy gadget mode */
 #define hydro_props_default_viscosity_alpha_max \
-  2.0f /* Values taken from (Price,2004), not used in legacy gadget mode */
+  2.0f /* values taken from (price,2004), not used in legacy gadget mode */
 #define hydro_props_default_viscosity_length \
   0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
+#endif /* ANARCHY_PU_SPH */
+
+/* Following values taken directly from the ANARCHY paper (Schaller+ 2015) */
+#define hydro_props_default_diffusion_alpha 0.0f
+#define hydro_props_default_diffusion_beta 0.01f
+#define hydro_props_default_diffusion_alpha_max 1.0f
+#define hydro_props_default_diffusion_alpha_min 0.0f
 
 /**
  * @brief Initialize the global properties of the hydro scheme.
@@ -143,6 +161,21 @@ void hydro_props_init(struct hydro_props *p,
   p->viscosity.length = parser_get_opt_param_float(
       params, "SPH:viscosity_length", hydro_props_default_viscosity_length);
 
+  /* Same for the thermal diffusion parameters */
+  p->diffusion.alpha = parser_get_opt_param_float(
+      params, "SPH:diffusion_alpha", hydro_props_default_diffusion_alpha);
+
+  p->diffusion.beta = parser_get_opt_param_float(
+      params, "SPH:diffusion_beta", hydro_props_default_diffusion_beta);
+
+  p->diffusion.alpha_max =
+      parser_get_opt_param_float(params, "SPH:diffusion_alpha_max",
+                                 hydro_props_default_diffusion_alpha_max);
+
+  p->diffusion.alpha_min =
+      parser_get_opt_param_float(params, "SPH:diffusion_alpha_min",
+                                 hydro_props_default_diffusion_alpha_min);
+
   /* Compute the initial energy (Note the temp. read is in internal units) */
   /* u_init = k_B T_init / (mu m_p (gamma - 1)) */
   double u_init = phys_const->const_boltzmann_k / phys_const->const_proton_mass;
@@ -312,6 +345,11 @@ void hydro_props_init_no_hydro(struct hydro_props *p) {
   p->viscosity.alpha_max = hydro_props_default_viscosity_alpha_max;
   p->viscosity.alpha_min = hydro_props_default_viscosity_alpha_min;
   p->viscosity.length = hydro_props_default_viscosity_length;
+
+  p->diffusion.alpha = hydro_props_default_diffusion_alpha;
+  p->diffusion.beta = hydro_props_default_diffusion_beta;
+  p->diffusion.alpha_max = hydro_props_default_diffusion_alpha_max;
+  p->diffusion.alpha_min = hydro_props_default_diffusion_alpha_min;
 }
 
 /**
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index 5ee6a22d2cf1c22f99e50a9254ef323d333d2a10..c63ecd9af777b6420858ee2a732d3012c123b209 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -105,6 +105,24 @@ struct hydro_props {
     /*! The decay length of the artificial viscosity (used in M&M, etc.) */
     float length;
   } viscosity;
+
+  /*! Thermal diffusion parameters */
+  struct {
+
+    /*! Initialisation value, or the case for constant thermal diffusion coeffs
+     */
+    float alpha;
+
+    /*! Tuning parameter for speed of ramp up/down */
+    float beta;
+
+    /*! Maximal value for alpha_diff */
+    float alpha_max;
+
+    /*! Minimal value for alpha_diff */
+    float alpha_min;
+
+  } diffusion;
 };
 
 void hydro_props_print(const struct hydro_props *p);
diff --git a/src/kernel_hydro.h b/src/kernel_hydro.h
index aac06c19ce39c647ba7211f85ac0a849365d126f..3d5ec6ac84a77941739f5d3b57ed0340c831c061 100644
--- a/src/kernel_hydro.h
+++ b/src/kernel_hydro.h
@@ -124,20 +124,27 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
 #define kernel_name "Wendland C2"
 #define kernel_degree 5 /* Degree of the polynomial */
 #define kernel_ivals 1  /* Number of branches */
+#if defined(HYDRO_DIMENSION_1D)
+/* Wendland C* have different form in 1D than 2D/3D */
+#define kernel_gamma ((float)(1.620185))
+#define kernel_constant ((float)(5. / 4.))
+static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
+    __attribute__((aligned(16))) = {
+        0.f, -3.f, 8.f, -6.f, 0.f, 1.f, /* 0 < u < 1 */
+        0.f, 0.f,  0.f, 0.f,  0.f, 0.f};
+#else
 #if defined(HYDRO_DIMENSION_3D)
 #define kernel_gamma ((float)(1.936492))
 #define kernel_constant ((float)(21. * M_1_PI / 2.))
 #elif defined(HYDRO_DIMENSION_2D)
 #define kernel_gamma ((float)(1.897367))
 #define kernel_constant ((float)(7. * M_1_PI))
-#elif defined(HYDRO_DIMENSION_1D)
-#error "Wendland C2 kernel not defined in 1D."
 #endif
 static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
     __attribute__((aligned(16))) = {
         4.f, -15.f, 20.f, -10.f, 0.f, 1.f,  /* 0 < u < 1 */
         0.f, 0.f,   0.f,  0.f,   0.f, 0.f}; /* 1 < u */
-
+#endif
 /* ------------------------------------------------------------------------- */
 #elif defined(WENDLAND_C4_KERNEL)
 
diff --git a/src/part.h b/src/part.h
index 64babf4a37696d7cb49b4804ee77773b4e1981fc..a0a09457039daa9ba59b7378638db12b5debafe2 100644
--- a/src/part.h
+++ b/src/part.h
@@ -75,6 +75,10 @@
 #elif defined(PLANETARY_SPH)
 #include "./hydro/Planetary/hydro_part.h"
 #define hydro_need_extra_init_loop 0
+#elif defined(ANARCHY_PU_SPH)
+#include "./hydro/AnarchyPU/hydro_part.h"
+#define hydro_need_extra_init_loop 0
+#define EXTRA_HYDRO_LOOP
 #else
 #error "Invalid choice of SPH variant"
 #endif
diff --git a/src/tools.c b/src/tools.c
index ca531671a7a2522eab760c1eb4896a6bd522a073..45b04e35bb4713a5b73158839e642d0796aa94a7 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -259,6 +259,85 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
   }
 }
 
+#ifdef EXTRA_HYDRO_LOOP
+void pairs_all_gradient(struct runner *r, struct cell *ci, struct cell *cj) {
+
+  float r2, hi, hj, hig2, hjg2, dx[3];
+  struct part *pi, *pj;
+  const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
+  const struct engine *e = r->e;
+  const struct cosmology *cosmo = e->cosmology;
+  const float a = cosmo->a;
+  const float H = cosmo->H;
+
+  /* Implements a double-for loop and checks every interaction */
+  for (int i = 0; i < ci->hydro.count; ++i) {
+
+    pi = &ci->hydro.parts[i];
+    hi = pi->h;
+    hig2 = hi * hi * kernel_gamma2;
+
+    /* Skip inactive particles. */
+    if (!part_is_active(pi, e)) continue;
+
+    for (int j = 0; j < cj->hydro.count; ++j) {
+
+      pj = &cj->hydro.parts[j];
+      hj = pj->h;
+      hjg2 = hj * hj * kernel_gamma2;
+
+      /* Pairwise distance */
+      r2 = 0.0f;
+      for (int k = 0; k < 3; k++) {
+        dx[k] = ci->hydro.parts[i].x[k] - cj->hydro.parts[j].x[k];
+        dx[k] = nearest(dx[k], dim[k]);
+        r2 += dx[k] * dx[k];
+      }
+
+      /* Hit or miss? */
+      if (r2 < hig2 && !part_is_inhibited(pj, e)) {
+
+        /* Interact */
+        runner_iact_nonsym_gradient(r2, dx, hi, hj, pi, pj, a, H);
+      }
+    }
+  }
+
+  /* Reverse double-for loop and checks every interaction */
+  for (int j = 0; j < cj->hydro.count; ++j) {
+
+    pj = &cj->hydro.parts[j];
+    hj = pj->h;
+    hjg2 = hj * hj * kernel_gamma2;
+
+    /* Skip inactive particles. */
+    if (!part_is_active(pj, e)) continue;
+
+    for (int i = 0; i < ci->hydro.count; ++i) {
+
+      pi = &ci->hydro.parts[i];
+      hi = pi->h;
+      hig2 = hi * hi * kernel_gamma2;
+
+      /* Pairwise distance */
+      r2 = 0.0f;
+      for (int k = 0; k < 3; k++) {
+        dx[k] = cj->hydro.parts[j].x[k] - ci->hydro.parts[i].x[k];
+        dx[k] = nearest(dx[k], dim[k]);
+        r2 += dx[k] * dx[k];
+      }
+
+      /* Hit or miss? */
+      if (r2 < hjg2 && !part_is_inhibited(pi, e)) {
+
+        /* Interact */
+        runner_iact_nonsym_gradient(r2, dx, hj, pi->h, pj, pi, a, H);
+      }
+    }
+  }
+}
+#endif /* EXTRA_HDYRO_LOOP */
+
 void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
 
   float r2, hi, hj, hig2, hjg2, dx[3];
@@ -460,6 +539,59 @@ void self_all_density(struct runner *r, struct cell *ci) {
   }
 }
 
+#ifdef EXTRA_HYDRO_LOOP
+void self_all_gradient(struct runner *r, struct cell *ci) {
+  float r2, hi, hj, hig2, hjg2, dxi[3];  //, dxj[3];
+  struct part *pi, *pj;
+  const struct engine *e = r->e;
+  const struct cosmology *cosmo = e->cosmology;
+  const float a = cosmo->a;
+  const float H = cosmo->H;
+
+  /* Implements a double-for loop and checks every interaction */
+  for (int i = 0; i < ci->hydro.count; ++i) {
+
+    pi = &ci->hydro.parts[i];
+    hi = pi->h;
+    hig2 = hi * hi * kernel_gamma2;
+
+    for (int j = i + 1; j < ci->hydro.count; ++j) {
+
+      pj = &ci->hydro.parts[j];
+      hj = pj->h;
+      hjg2 = hj * hj * kernel_gamma2;
+
+      if (pi == pj) continue;
+
+      /* Pairwise distance */
+      r2 = 0.0f;
+      for (int k = 0; k < 3; k++) {
+        dxi[k] = ci->hydro.parts[i].x[k] - ci->hydro.parts[j].x[k];
+        r2 += dxi[k] * dxi[k];
+      }
+
+      /* Hit or miss? */
+      if (r2 < hig2 && part_is_active(pi, e) && !part_is_inhibited(pj, e)) {
+
+        /* Interact */
+        runner_iact_nonsym_gradient(r2, dxi, hi, hj, pi, pj, a, H);
+      }
+
+      /* Hit or miss? */
+      if (r2 < hjg2 && part_is_active(pj, e) && !part_is_inhibited(pi, e)) {
+
+        dxi[0] = -dxi[0];
+        dxi[1] = -dxi[1];
+        dxi[2] = -dxi[2];
+
+        /* Interact */
+        runner_iact_nonsym_gradient(r2, dxi, hj, hi, pj, pi, a, H);
+      }
+    }
+  }
+}
+#endif /* EXTRA_HYDRO_LOOP */
+
 void self_all_force(struct runner *r, struct cell *ci) {
   float r2, hi, hj, hig2, hjg2, dxi[3];  //, dxj[3];
   struct part *pi, *pj;
diff --git a/src/tools.h b/src/tools.h
index 22eba1519ebf80673cb2d8540791e6b4d092bab0..a34904bcbb7af1ce15408a376f957f0f72cd327c 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -38,6 +38,8 @@ void pairs_single_density(double *dim, long long int pid,
 
 void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj);
 void self_all_density(struct runner *r, struct cell *ci);
+void pairs_all_gradient(struct runner *r, struct cell *ci, struct cell *cj);
+void self_all_gradient(struct runner *r, struct cell *ci);
 void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj);
 void self_all_force(struct runner *r, struct cell *ci);
 void pairs_all_stars_density(struct runner *r, struct cell *ci,
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 5b518970ea118c98a8354e816f86ecd16a5f85cf..7a05be68fd6aeb0ac2bb13d28209c9022e97b633 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -112,7 +112,7 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
 #elif defined(DEFAULT_SPH)
   part->u = pressure / (hydro_gamma_minus_one * density);
 #elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
-    defined(HOPKINS_PU_SPH_MONAGHAN)
+    defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH)
   part->u = pressure / (hydro_gamma_minus_one * density);
 #elif defined(PLANETARY_SPH)
   part->u = pressure / (hydro_gamma_minus_one * density);
@@ -401,10 +401,13 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
             main_cell->hydro.parts[pid].v[0], main_cell->hydro.parts[pid].v[1],
             main_cell->hydro.parts[pid].v[2], main_cell->hydro.parts[pid].h,
             hydro_get_comoving_density(&main_cell->hydro.parts[pid]),
-#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) ||   \
-    defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH) || \
-    defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN)
+#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) ||              \
+    defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH) ||            \
+    defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN) || \
+    defined(ANARCHY_PU_SPH)
             0.f,
+#elif defined(ANARCHY_PU_SPH)
+            main_cell->hydro.parts[pid].viscosity.div_v,
 #else
             main_cell->hydro.parts[pid].density.div_v,
 #endif
@@ -427,6 +430,9 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
     defined(HOPKINS_PU_SPH_MONAGHAN)
             main_cell->hydro.parts[pid].force.v_sig, 0.f,
             main_cell->hydro.parts[pid].u_dt
+#elif defined(ANARCHY_PU_SPH)
+            main_cell->hydro.parts[pid].viscosity.v_sig, 0.f,
+            main_cell->hydro.parts[pid].u_dt
 #else
             0.f, 0.f, 0.f
 #endif
@@ -460,9 +466,16 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
 void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
                                    struct cell *cj);
 void runner_doself1_branch_density(struct runner *r, struct cell *ci);
+#ifdef EXTRA_HYDRO_LOOP
+void runner_dopair1_branch_gradient(struct runner *r, struct cell *ci,
+                                    struct cell *cj);
+void runner_doself1_branch_gradient(struct runner *r, struct cell *ci);
+#endif /* EXTRA_HYDRO LOOP */
 void runner_dopair2_branch_force(struct runner *r, struct cell *ci,
                                  struct cell *cj);
 void runner_doself2_branch_force(struct runner *r, struct cell *ci);
+void runner_doself2_force(struct runner *r, struct cell *ci);
+void runner_doself2_force_vec(struct runner *r, struct cell *ci);
 
 /* And go... */
 int main(int argc, char *argv[]) {
@@ -588,7 +601,7 @@ int main(int argc, char *argv[]) {
   hp.eta_neighbours = h;
   hp.h_tolerance = 1e0;
   hp.h_max = FLT_MAX;
-  hp.max_smoothing_iterations = 1;
+  hp.max_smoothing_iterations = 10;
   hp.CFL_condition = 0.1;
 
   struct engine engine;
@@ -674,7 +687,7 @@ int main(int argc, char *argv[]) {
     cache_init(&runner.cj_cache, 512);
 #endif
 
-    /* Run all the pairs (only once !)*/
+    /* Run all the  (only once !)*/
     for (int i = 0; i < 5; i++) {
       for (int j = 0; j < 5; j++) {
         for (int k = 0; k < 5; k++) {
@@ -711,6 +724,51 @@ int main(int argc, char *argv[]) {
     /* Ghost to finish everything on the central cells */
     for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j], 0);
 
+#ifdef EXTRA_HYDRO_LOOP
+    /* We need to do the gradient loop and the extra ghost! */
+    message(
+        "Extra hydro loop detected, running gradient loop in test125cells.");
+
+    /* Run all the pairs (only once !)*/
+    for (int i = 0; i < 5; i++) {
+      for (int j = 0; j < 5; j++) {
+        for (int k = 0; k < 5; k++) {
+
+          struct cell *ci = cells[i * 25 + j * 5 + k];
+
+          for (int ii = -1; ii < 2; ii++) {
+            int iii = i + ii;
+            if (iii < 0 || iii >= 5) continue;
+            iii = (iii + 5) % 5;
+            for (int jj = -1; jj < 2; jj++) {
+              int jjj = j + jj;
+              if (jjj < 0 || jjj >= 5) continue;
+              jjj = (jjj + 5) % 5;
+              for (int kk = -1; kk < 2; kk++) {
+                int kkk = k + kk;
+                if (kkk < 0 || kkk >= 5) continue;
+                kkk = (kkk + 5) % 5;
+
+                struct cell *cj = cells[iii * 25 + jjj * 5 + kkk];
+
+                if (cj > ci) runner_dopair1_branch_gradient(&runner, ci, cj);
+              }
+            }
+          }
+        }
+      }
+    }
+
+    /* And now the self-interaction for the central cells */
+    for (int j = 0; j < 27; ++j)
+      runner_doself1_branch_gradient(&runner, inner_cells[j]);
+
+    /* Extra ghost to finish everything on the central cells */
+    for (int j = 0; j < 27; ++j)
+      runner_do_extra_ghost(&runner, inner_cells[j], 0);
+
+#endif /* EXTRA_HYDRO_LOOP */
+
       /* Do the force calculation */
 
 #ifdef WITH_VECTORIZATION
@@ -779,11 +837,11 @@ int main(int argc, char *argv[]) {
 
   ticks self_time = timings[26];
 
-  message("Corner calculations took       : %15lli ticks.", corner_time / runs);
-  message("Edge calculations took         : %15lli ticks.", edge_time / runs);
-  message("Face calculations took         : %15lli ticks.", face_time / runs);
-  message("Self calculations took         : %15lli ticks.", self_time / runs);
-  message("SWIFT calculation took         : %15lli ticks.", time / runs);
+  message("Corner calculations took:     %15lli ticks.", corner_time / runs);
+  message("Edge calculations took:       %15lli ticks.", edge_time / runs);
+  message("Face calculations took:       %15lli ticks.", face_time / runs);
+  message("Self calculations took:       %15lli ticks.", self_time / runs);
+  message("SWIFT calculation took:       %15lli ticks.", time / runs);
 
   for (int j = 0; j < 125; ++j)
     reset_particles(cells[j], &space.hs, vel, press, size, rho);
@@ -840,6 +898,48 @@ int main(int argc, char *argv[]) {
   /* Ghost to finish everything on the central cells */
   for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j], 0);
 
+#ifdef EXTRA_HYDRO_LOOP
+  /* We need to do the gradient loop and the extra ghost! */
+
+  /* Run all the pairs (only once !)*/
+  for (int i = 0; i < 5; i++) {
+    for (int j = 0; j < 5; j++) {
+      for (int k = 0; k < 5; k++) {
+
+        struct cell *ci = cells[i * 25 + j * 5 + k];
+
+        for (int ii = -1; ii < 2; ii++) {
+          int iii = i + ii;
+          if (iii < 0 || iii >= 5) continue;
+          iii = (iii + 5) % 5;
+          for (int jj = -1; jj < 2; jj++) {
+            int jjj = j + jj;
+            if (jjj < 0 || jjj >= 5) continue;
+            jjj = (jjj + 5) % 5;
+            for (int kk = -1; kk < 2; kk++) {
+              int kkk = k + kk;
+              if (kkk < 0 || kkk >= 5) continue;
+              kkk = (kkk + 5) % 5;
+
+              struct cell *cj = cells[iii * 25 + jjj * 5 + kkk];
+
+              if (cj > ci) pairs_all_gradient(&runner, ci, cj);
+            }
+          }
+        }
+      }
+    }
+  }
+
+  /* And now the self-interaction for the central cells */
+  for (int j = 0; j < 27; ++j) self_all_gradient(&runner, inner_cells[j]);
+
+  /* Extra ghost to finish everything on the central cells */
+  for (int j = 0; j < 27; ++j)
+    runner_do_extra_ghost(&runner, inner_cells[j], 0);
+
+#endif /* EXTRA_HYDRO_LOOP */
+
   /* Do the force calculation */
 
   /* Do the pairs (for the central 27 cells) */
@@ -864,7 +964,7 @@ int main(int argc, char *argv[]) {
   const ticks toc = getticks();
 
   /* Output timing */
-  message("Brute force calculation took : %15lli ticks.", toc - tic);
+  message("Brute force calculation took: %15lli ticks.", toc - tic);
 
   sprintf(outputFileName, "brute_force_125_%.150s.dat",
           outputFileNameExtension);
diff --git a/tests/test27cells.c b/tests/test27cells.c
index d100e2e30f0bb1452b3366ddde51dbae0575d67a..cc34f503304feb56799a2d31baa3416b940202d3 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -275,7 +275,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
             hydro_get_comoving_density(&main_cell->hydro.parts[pid]),
 #if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
             0.f,
-#elif defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN)
+#elif defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN) || \
+    defined(ANARCHY_PU_SPH)
             main_cell->hydro.parts[pid].density.pressure_bar_dh,
 #else
             main_cell->hydro.parts[pid].density.rho_dh,
@@ -288,6 +289,12 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
             main_cell->hydro.parts[pid].density.rot_v[0],
             main_cell->hydro.parts[pid].density.rot_v[1],
             main_cell->hydro.parts[pid].density.rot_v[2]
+#elif defined(ANARCHY_PU_SPH)
+            /* this is required because of the variable AV scheme */
+            main_cell->hydro.parts[pid].viscosity.div_v,
+            main_cell->hydro.parts[pid].density.rot_v[0],
+            main_cell->hydro.parts[pid].density.rot_v[1],
+            main_cell->hydro.parts[pid].density.rot_v[2]
 #else
             0., 0., 0., 0.
 #endif
@@ -327,6 +334,12 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
               cj->hydro.parts[pjd].density.rot_v[0],
               cj->hydro.parts[pjd].density.rot_v[1],
               cj->hydro.parts[pjd].density.rot_v[2]
+#elif defined(ANARCHY_PU_SPH)
+              /* this is required because of the variable AV scheme */
+              cj->hydro.parts[pjd].viscosity.div_v,
+              cj->hydro.parts[pjd].density.rot_v[0],
+              cj->hydro.parts[pjd].density.rot_v[1],
+              cj->hydro.parts[pjd].density.rot_v[2]
 #else
               0., 0., 0., 0.
 #endif
diff --git a/tests/testActivePair.c b/tests/testActivePair.c
index 402a9a7ac416a0b2651f628eb32988d8ad62a14f..54a3189b89d9de757bf340bf759db5b40f947174 100644
--- a/tests/testActivePair.c
+++ b/tests/testActivePair.c
@@ -112,7 +112,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
 #if defined(GADGET2_SPH)
         part->entropy = 1.f;
 #elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
-    defined(HOPKINS_PU_SPH_MONAGHAN)
+    defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH)
         part->u = 1.f;
 #elif defined(HOPKINS_PE_SPH)
         part->entropy = 1.f;
@@ -212,7 +212,8 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo,
     p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
     p->density.wcount_dh = 0.f;
 #endif /* PRESSURE-ENTROPY */
-#if defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN)
+#if defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN) || \
+    defined(ANARCHY_PU_SPH)
     p->rho = 1.f;
     p->pressure_bar = 0.6666666;
     p->density.rho_dh = 0.f;
@@ -220,6 +221,13 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo,
     p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
     p->density.wcount_dh = 0.f;
 #endif /* PRESSURE-ENERGY */
+#if defined(ANARCHY_PU_SPH)
+    /* Initialise viscosity variables */
+    p->viscosity.alpha = 0.8;
+    p->viscosity.div_v = 0.f;
+    p->viscosity.div_v_previous_step = 0.f;
+    p->viscosity.v_sig = hydro_get_comoving_soundspeed(p);
+#endif /* ANARCHY_PU_SPH viscosity variables */
 
     /* And prepare for a round of force tasks. */
     hydro_prepare_force(p, xp, cosmo, hydro_props, 0.);
diff --git a/tests/testCbrt.c b/tests/testCbrt.c
index 3663e0e19ad2a5ad35d67703e00f5c0309a3eb00..bba379902b2bbc16bd49a5bbba0917100b4d60a7 100644
--- a/tests/testCbrt.c
+++ b/tests/testCbrt.c
@@ -38,7 +38,9 @@ int main(int argc, char *argv[]) {
   clocks_set_cpufreq(cpufreq);
 
   /* Choke on FP-exceptions */
+#ifdef HAVE_FE_ENABLE_EXCEPT
   feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+#endif
 
   /* Some constants for this test. */
   const int num_vals = 200000000;
@@ -59,8 +61,8 @@ int main(int argc, char *argv[]) {
   for (int k = 0; k < num_vals; k++) {
     const double exact = cbrt(data[k]);  // computed in double just to be sure.
     const float ours = 1.0f / icbrtf(data[k]);
-    const float err_abs = fabsf(exact - ours);
-    const float err_rel = 0.5f * fabsf(exact - ours) / (exact + ours);
+    const float err_abs = fabs(exact - ours);
+    const float err_rel = 0.5f * fabs(exact - ours) / (exact + ours);
 
     if (err_rel > err_rel_tol && data[k] != 0.f)
       error(
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index dae55a337642e1616e94119263ff8f1c2a617c89..e14fddd640764c7e22a217fb483791494ba4fae0 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -111,9 +111,10 @@ struct part *make_particles(size_t count, double *offset, double spacing,
  */
 void prepare_force(struct part *parts, size_t count) {
 
-#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH) && \
-    !defined(MINIMAL_SPH) && !defined(PLANETARY_SPH) &&   \
-    !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN)
+#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH) &&            \
+    !defined(MINIMAL_SPH) && !defined(PLANETARY_SPH) &&              \
+    !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN) && \
+    !defined(ANARCHY_PU_SPH)
   struct part *p;
   for (size_t i = 0; i < count; ++i) {
     p = &parts[i];
@@ -154,7 +155,7 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
 #elif defined(DEFAULT_SPH)
           p->force.v_sig, 0.f, p->force.u_dt
 #elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
-    defined(HOPKINS_PU_SPH_MONAGHAN)
+    defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH)
           p->force.v_sig, 0.f, p->u_dt
 #else
           0.f, 0.f, 0.f
@@ -558,7 +559,8 @@ void test_force_interactions(struct part test_part, struct part *parts,
       vizq[i] = pi_vec.v[2];
       rhoiq[i] = pi_vec.rho;
       grad_hiq[i] = pi_vec.force.f;
-#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN)
+#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN) && \
+    !defined(ANARCHY_PU_SPH)
       pOrhoi2q[i] = pi_vec.force.P_over_rho2;
 #endif
       balsaraiq[i] = pi_vec.force.balsara;
@@ -571,7 +573,8 @@ void test_force_interactions(struct part test_part, struct part *parts,
       vjzq[i] = pj_vec[i].v[2];
       rhojq[i] = pj_vec[i].rho;
       grad_hjq[i] = pj_vec[i].force.f;
-#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN)
+#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN) && \
+    !defined(ANARCHY_PU_SPH)
       pOrhoj2q[i] = pj_vec[i].force.P_over_rho2;
 #endif
       balsarajq[i] = pj_vec[i].force.balsara;
@@ -653,7 +656,8 @@ void test_force_interactions(struct part test_part, struct part *parts,
     VEC_HADD(a_hydro_zSum, piq[0]->a_hydro[2]);
     VEC_HADD(h_dtSum, piq[0]->force.h_dt);
     VEC_HMAX(v_sigSum, piq[0]->force.v_sig);
-#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN)
+#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN) && \
+    !defined(ANARCHY_PU_SPH)
     VEC_HADD(entropy_dtSum, piq[0]->entropy_dt);
 #endif
 
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
index ce1e2e9354c4d59a6e58619d43b743864ed38585..1f0849bb9093948fa68d88984c285c44b403ba79 100644
--- a/tests/testSymmetry.c
+++ b/tests/testSymmetry.c
@@ -177,6 +177,38 @@ void test(void) {
     error("Particles 'pj' do not match after density (byte = %d)", j_not_ok);
   }
 
+    /* --- Test the gradient loop --- */
+#ifdef EXTRA_HYDRO_LOOP
+
+  /* Call the symmetric version */
+  runner_iact_gradient(r2, dx, pi.h, pj.h, &pi, &pj, a, H);
+
+  /* Call the non-symmetric version */
+  runner_iact_nonsym_gradient(r2, dx, pi2.h, pj2.h, &pi2, &pj2, a, H);
+  dx[0] = -dx[0];
+  dx[1] = -dx[1];
+  dx[2] = -dx[2];
+  runner_iact_nonsym_gradient(r2, dx, pj2.h, pi2.h, &pj2, &pi2, a, H);
+
+  i_not_ok = memcmp((char *)&pi, (char *)&pi2, sizeof(struct part));
+  j_not_ok = memcmp((char *)&pj, (char *)&pj2, sizeof(struct part));
+
+  if (i_not_ok) {
+    printParticle_single(&pi, &xpi);
+    printParticle_single(&pi2, &xpi);
+    print_bytes(&pi, sizeof(struct part));
+    print_bytes(&pi2, sizeof(struct part));
+    error("Particles 'pi' do not match after gradient (byte = %d)", i_not_ok);
+  }
+  if (j_not_ok) {
+    printParticle_single(&pj, &xpj);
+    printParticle_single(&pj2, &xpj);
+    print_bytes(&pj, sizeof(struct part));
+    print_bytes(&pj2, sizeof(struct part));
+    error("Particles 'pj' do not match after gradient (byte = %d)", j_not_ok);
+  }
+#endif
+
   /* --- Test the force loop --- */
 
   /* Call the symmetric version */
diff --git a/theory/SPH/Flavours/anarchy.tex b/theory/SPH/Flavours/anarchy.tex
new file mode 100644
index 0000000000000000000000000000000000000000..5924f9438f9b553298b0d45a8e4d7ddae9167270
--- /dev/null
+++ b/theory/SPH/Flavours/anarchy.tex
@@ -0,0 +1,123 @@
+\section{ANARCHY-SPH}
+\label{sec:sph:anarchy}
+
+This section is loosely based on Dalla Vecchia (\textit{in prep.}), also described in section 2.2.2 of
+\cite{Schaller2015}.\\
+
+The version of ANARCHY that is currently implemented in \swift{} is based on the Pressure-Energy
+(P-U) SPH scheme, rather than the original Pressure-Entropy. This was chosen to make the
+implementation of sub-grid physics easier, as well as injection of energy more accurate as
+no iteration of the entropy-energy relation is required.
+
+ANARCH SPH comprises of:
+\begin{itemize}
+	\item Pressure-Energy SPH
+	\item \citet[][, henceforth C\&D]{cullen2010} variable artificial viscosity
+	\item A basic thermal diffusion term
+	\item The time-step limiter from \citet{durier2012}.
+\end{itemize}
+
+\subsection{Equations of Motion}
+
+The following smoothed quantities are required, and are calculated in the density loop:
+\begin{itemize}
+	\item $\rho_i = [h_i^{-n_d}]\sum_j m_j w_{i}$
+	\item $(d\rho/dh)_i = - [h_i^{-n_d - 1}]\sum_j m_j ( n_d * w_i + x_i \nabla_i w_i)$
+	\item $\bar{P}_i = [(\gamma - 1)h_i^{-n_d}]\sum_j m_j u_j w_{i}$
+	\item $(d\bar{P}/dh)_i = - [(\gamma - 1)h_i^{-n_d - 1}]\sum_j m_j u_j ( n_d * w_i + x_i \nabla_i w_i)$
+	\item $n_i = [h_i^{-n_d}]\sum_j w_{i}$
+	\item $(dn/dh)_i = - [h_i^{-n_d - 1}]\sum_j ( n_d * w_i + x_i \nabla_i w_i)$
+	\item $(\nabla \cdot \mathbf{v})_i = - [a^{-2} \rho_i^{-1} h^{-n_d - 1}]\sum_j m_j \mathbf{v}_{ij} \cdot \tilde{\mathbf{x}}_{ij} \nabla_i w_i$
+	% Think the cosmo factor entering here is wrong...
+	\item $(\nabla \times \mathbf{v})_i = - [a^{-2} \rho_i^{-1} h^{-n_d - 1} + Hn_d]\sum_j m_j \mathbf{v}_{ij} \times \tilde{\mathbf{x}}_{ij} \nabla_i w_i$
+\end{itemize}
+with quantities in square brackets added in {\tt hydro\_end\_density} and:
+\begin{itemize}
+	\item $h_i$ the smoothing length of particle $i$
+	\item $n_d$ the number of hydro dimensions
+	\item $m_j$ the mass of particle $j$
+	\item $w_{i}$ the dimensionless kernel evalulated at $x_i = r / h_i$
+	\item $r$ the interparticle separation of particles $i$ and $j$
+	\item $u_i$ the internal energy of the particle
+	\item $\gamma$ the ratio of specific heats
+	\item $\mathbf{v}_{ij}$ the difference between the velocities of particle $i$ and $j$
+	\item $\tilde{\mathbf{x}}_{ij}$ the unit vector connecting particles $i$ and $j$
+	\item $a$ the current scale factor
+	\item $H$ the current Hubble constant
+\end{itemize}
+
+The ANARCHY scheme requires a gradient loop, as intermediate smoothed quantities are required
+for the artificial viscosity and diffusion schemes. The following quatntities are calculated:
+\begin{itemize}
+	\item $v_{{\rm sig}, i} = \rm{max}(v_{\rm sig}, c_i + c_j - 3\mathbf{v}_{ij} \cdot \tilde{\mathbf{x}}_{ij})$
+	\item $\nabla^2 u_i = [2]\sum_j m_j \frac{u_i - u_j}{\rho_j} \frac{\nabla_i W_i}{r_{ij}}$
+\end{itemize}
+with quantities in square brackets added in {\tt hydro\_end\_gradient} and:
+\begin{itemize}
+	\item $v_{\rm sig}$ the siginal velocity	
+	\item $c_i$ the sound speed of particle $i$
+\end{itemize}
+
+In {\tt hydro\_prepare\_force}, the differential equations for the viscosity and
+diffusion schemes are integrated as follows. This includes some logic, so it is
+split into viscosity:
+\begin{itemize}
+	\item $\tau_i = h_i / (2 v_{{\rm sig}, i} \ell$
+	\item $\dot{\nabla \cdot \mathbf{v}_i} =
+	       \left({\nabla \cdot \mathbf{v}_i}(t+dt) - {\nabla \cdot \mathbf{v}_i}(t)\right)
+	       / dt$
+	\item $S_i = h_i^2 {\rm max}(0, -\dot{\nabla \cdot \mathbf{v}_i})$
+	\item $\alpha_{{\rm loc}, i} = \alpha_{\rm max} S_i / (S_i + v_{{\rm sig}, i}^2)$.
+\end{itemize}
+and diffusion:
+\begin{itemize}
+	\item $\dot{\tilde{\alpha}}_i = \beta h_i \frac{\nabla^2 u_i}{\sqrt{u_i}}$
+\end{itemize}
+where:
+\begin{itemize}
+	\item $\alpha_i$ is the viscosity coefficient
+	\item $\tilde{\alpha}_i$ is the diffusion coefficient
+	\item $\tau_i$ is the timescale for decay
+	\item $\ell$ is the viscosity length coefficient
+	\item $\beta$ is the diffusion length coefficient
+\end{itemize}
+The equations are then integrated as follows for viscosity:
+\begin{enumerate}
+	\item If $\alpha_{\rm loc} > \alpha_i$, update $\alpha_i$ to $\alpha_{\rm loc}$
+	      immediately.
+	\item Otherwise, decay the viscosity, with $\dot{\alpha}_i = (\alpha_{\rm loc} - \alpha_i) / \tau_i$.
+	      This equation is integrated with the same time-step as the velocity divergence derivative
+	      uses, and is the same time-step used for the cooling.
+	\item Finally, if $\alpha_i < \alpha_{\rm min}$, $\alpha_i$ is reset to that minimal
+	      value.
+\end{enumerate}
+and for diffusion:
+\begin{enumerate}
+	\item First, find the new diffusion coefficient, $\tilde{\alpha}_i(t+dt) = 
+	      \tilde{\alpha}_i(t) + \dot{\tilde{\alpha}}_i \cdot dt$, using the
+	      same time-step as for the viscosity.
+	\item If this is outside of the bounds set for the coefficient, set it
+	      to the respective upper or lower bound.
+\end{enumerate}
+The final force loop calculates the equations of motion for the particles ready for
+their time-integration. The following quantities are calculated:
+\begin{itemize}
+	\item $\mathbf{a}_{\rm hydro} = -\sum_j m_j u_i u_j (\gamma - 1)^2 \left(
+	       \frac{f_{ij}}{\bar{P}_i} \nabla_i W_i + \frac{f_{ji}}{\bar{P}_j} \nabla_j W_j\right)$
+	\item $\mathbf{a}_{\rm visc} = - \frac{1}{8}\sum_j (\alpha_i + \alpha_j) v_{{\rm sig}, i}
+	       \mu_{ij} (b_i + b_j) (\nabla_i W_i + \nabla_j W_j)/ (\rho_i + \rho_j)$
+	\item $\dot{u}_{ij, {\rm hydro}} = \sum_j m_j u_i u_j (\gamma - 1)^2
+	       \frac{f_{ij}}{\bar{P}_i} \nabla_i W_i$
+	\item $\dot{u}_{ij, {\rm visc}} = \frac{1}{2} \a_{\rm visc} (\mathbf{v}_{ij} \cdot \tilde{\mathbf{x}}_{ij} + r^2a^2 H)$
+	\item $v_{{\rm diff}, i} = {\rm max}(0, c_i + c_j + \mathbf{v}_{ij} \cdot \tilde{\mathbf{x}}_{ij} + r^2a^2 H)$
+	\item $\dot{u}_{ij, {\rm diff}} = \frac{1}{2}(\tilde{\alpha}_i + \tilde{\alpha}_j) a^{(3\gamma - 5)/2)}
+	       v_{{\rm diff}, i} (u_i - u_j) (\nabla_i W_i + \nabla_j W_j)/ (\rho_i + \rho_j) $
+	\item $\dot{u}_i = \sum_j \dot{u}_{ij, {\rm hydro}} +  \dot{u}_{ij, {\rm visc}} + \dot{u}_{ij, {\rm diff}}$
+\end{itemize}
+where:
+\begin{itemize}
+	\item $f_{ij}$ are the variable smoothing length correction factors
+	\item $b_i$ is the Balsara switch for particle $i$
+	\item $\mu_{ij} = a^{(3\gamma - 5)/2) {\rm min}(\mathbf{v}_{ij} \cdot \tilde{\mathbf{x}}_{ij} + r^2a^2 H, 0)$
+\end{itemize}
+
diff --git a/theory/SPH/Flavours/sph_flavours.tex b/theory/SPH/Flavours/sph_flavours.tex
index 81ac2153ed23f29f14345e0377774548420c84c9..d84bdcc0b42129e9d5008051dd6e0e212e4e9463 100644
--- a/theory/SPH/Flavours/sph_flavours.tex
+++ b/theory/SPH/Flavours/sph_flavours.tex
@@ -626,8 +626,3 @@ The rest of the artificial viscosity implementation, including the
 \citet{Balsara1995} switch, is the same - just with $\alpha \rightarrow
 \alpha_{ij}$.
 
-\subsection{Anarchy SPH}
-Dalla Vecchia (\textit{in prep.}), also described in section 2.2.2 of
-\cite{Schaller2015}.\\
-\label{sec:sph:anarchy}
-\tbd 
diff --git a/theory/SPH/Flavours/sph_flavours_standalone.tex b/theory/SPH/Flavours/sph_flavours_standalone.tex
index 20c9f1451c2499d661bfbb1022bfd34c02fde4dd..7cc92fdb438ea09916dfbd12bfb554ddb9fc2ecd 100644
--- a/theory/SPH/Flavours/sph_flavours_standalone.tex
+++ b/theory/SPH/Flavours/sph_flavours_standalone.tex
@@ -24,6 +24,7 @@
 
 \maketitle
 \input{sph_flavours}
+\input{anarchy}
 
 \bibliographystyle{mnras}
 \bibliography{./bibliography}
diff --git a/theory/SPH/swift_sph.tex b/theory/SPH/swift_sph.tex
index e9c185c3cd0b845bff75be2092846bffbdcfd1a9..51ab6c3e49ae6b74d8b63590caeadd962cd6e4d5 100644
--- a/theory/SPH/swift_sph.tex
+++ b/theory/SPH/swift_sph.tex
@@ -34,6 +34,7 @@
 
 \section{SPH flavours}
 \input{Flavours/sph_flavours}
+\input{Flavours/anarchy}
 
 \bibliographystyle{mnras}
 \bibliography{./bibliography}