diff --git a/doc/RTD/source/CommandLineOptions/index.rst b/doc/RTD/source/CommandLineOptions/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..454a2e13433bd34349eea1337fdde854e7c1a4c2
--- /dev/null
+++ b/doc/RTD/source/CommandLineOptions/index.rst
@@ -0,0 +1,6 @@
+.. Command line options
+   Matthieu Schaller, 21st October 2018
+
+Command line options
+====================
+
diff --git a/doc/RTD/source/GettingStarted/index.rst b/doc/RTD/source/GettingStarted/index.rst
index 36de8ea740490c16bc9d6b69d871290e80dc2091..2086bcfb4af0ac1b7bbc24c34caa85fa1ebec498 100644
--- a/doc/RTD/source/GettingStarted/index.rst
+++ b/doc/RTD/source/GettingStarted/index.rst
@@ -20,7 +20,6 @@ and keep on your desk.
    running_example
    runtime_options
    configuration_options
-   parameter_file
    what_about_mpi
    running_on_large_systems
    special_modes
diff --git a/doc/RTD/source/HydroSchemes/hopkins_sph.rst b/doc/RTD/source/HydroSchemes/hopkins_sph.rst
index bcc51e0ad96b18956f1c8e54f7bf2bf3b352c138..e4f1479230df96eabaa1fe16994960059858613b 100644
--- a/doc/RTD/source/HydroSchemes/hopkins_sph.rst
+++ b/doc/RTD/source/HydroSchemes/hopkins_sph.rst
@@ -28,3 +28,9 @@ scheme it includes a Monaghan AV scheme and a Balsara switch.
 .. code-block:: bash
    
    ./configure --with-hydro="pressure-energy"
+
+Both of the above schemes use a very simple, fixed artificial viscosity, only
+the ``SPH:viscosity_alpha`` parameter has any effect for this scheme. This will
+change the strength of the artificial viscosity throughout the simulation, and
+has a default of 0.8.
+
diff --git a/doc/RTD/source/HydroSchemes/minimal_sph.rst b/doc/RTD/source/HydroSchemes/minimal_sph.rst
index 1a16a23360aaba8b28920150af0d4f4b05c74c2f..bbcbe026b56381c007f58920f31115f9f9160d05 100644
--- a/doc/RTD/source/HydroSchemes/minimal_sph.rst
+++ b/doc/RTD/source/HydroSchemes/minimal_sph.rst
@@ -10,11 +10,17 @@ Minimal (Density-Energy) SPH
    :caption: Contents:
 
 This scheme is a textbook implementation of Density-Energy SPH, and can be used
-as a pedagogical example. It also implements a Monaghan AV scheme, like the
-GADGET-2 scheme. It uses very similar equations, but differs in implementation
-details; namely it tracks the internal energy \(u\) as the thermodynamic
-variable, rather than entropy \(A\). To use the minimal scheme, use
+as a pedagogical example. It also implements a Monaghan AV scheme with a
+Balsara switch, like the GADGET-2 scheme. It uses very similar equations, but
+differs in implementation details; namely it tracks the internal energy \(u\)
+as the thermodynamic variable, rather than entropy \(A\). To use the minimal
+scheme, use
 
 .. code-block:: bash
 
     ./configure --with-hydro="minimal"
+
+As it uses a very simple, fixed artificial viscosity, only the
+``SPH:viscosity_alpha`` parameter has any effect for this scheme. This will
+change the strength of the artificial viscosity throughout the simulation,
+and has a default of 0.8.
diff --git a/doc/RTD/source/HydroSchemes/traditional_sph.rst b/doc/RTD/source/HydroSchemes/traditional_sph.rst
index c69ea5f60644119b8590414ffe00a75246de49a6..455e8bebe516bd9be9f6df889f1ead2088ca94d2 100644
--- a/doc/RTD/source/HydroSchemes/traditional_sph.rst
+++ b/doc/RTD/source/HydroSchemes/traditional_sph.rst
@@ -15,3 +15,8 @@ a Monaghan artificial viscosity scheme and Balsara switch.
 To use this hydro scheme, you need no extra configuration options -- it is the
 default!
 
+As it uses a very simple, fixed artificial viscosity, only the
+``SPH:viscosity_alpha`` parameter has any effect for this scheme. This will
+change the strength of the artificial viscosity throughout the simulation,
+and has a default of 0.8.
+
diff --git a/doc/RTD/source/InitialConditions/index.rst b/doc/RTD/source/InitialConditions/index.rst
index 19ead066c4efb799e6587884a952baf488f794bf..f9de5364ca49bfb602a78e875ea5dd33633a1321 100644
--- a/doc/RTD/source/InitialConditions/index.rst
+++ b/doc/RTD/source/InitialConditions/index.rst
@@ -176,7 +176,9 @@ as peculiar velocities divided by ``sqrt(a)``. This can be undone by swicthing
 on the parameter ``InitialConditions:cleanup_velocity_factors`` in the
 :ref:`Parameter_File_label`.
 
-     
+
+.. _ICs_units_label:
+
 Optional Components
 -------------------
 
diff --git a/doc/RTD/source/ParameterFiles/index.rst b/doc/RTD/source/ParameterFiles/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..3d865dac22f3a2083a3627c2ae3813dec0490ff5
--- /dev/null
+++ b/doc/RTD/source/ParameterFiles/index.rst
@@ -0,0 +1,237 @@
+.. Parameter Files
+   Matthieu Schaller, 21st October 2018
+
+Parameter Files
+===============
+
+File format and basic information
+---------------------------------
+
+The parameter file uses a format similar to the `YAML format
+<https://en.wikipedia.org/wiki/YAML>`_ but reduced to only the
+elements required for the SWIFT parameters. Options are given by a
+name followed by a column and the value of the parameter:
+
+.. code:: YAML
+
+   ICs:        santa_barbara.hdf5	  
+   dt_max:     1.5
+   shift:      [2., 4., 5.]
+
+Comments can be inserted anywhere and start with a hash:
+
+.. code:: YAML
+
+   # Descrption of the physics
+   viscosity_alpha:     2.0
+   dt_max:              1.5     # seconds
+
+A typical SWIFT parameter file is split into multiple sections that
+may or may not be present depending on the different configuration
+options. The sections start with a label and can contain any number of
+parameters:
+
+.. code:: YAML
+
+   Cosmology:    # Planck13
+     Omega_m:        0.307
+     Omega_lambda:   0.693
+     Omega_b:        0.0455
+     h:              0.6777
+     a_begin:        0.0078125     # z = 127
+
+The options can be integer values, floating point numbers, characters
+or strings. If SWIFT expects a number and string is given, an error
+will be raised. The code can also read an array of values:
+
+.. code:: YAML
+
+   shift:  [2., 4., 5.]
+	  
+Some options in the parameter file are optional and
+when not provided, SWIFT will run with the default value. However, if
+a compulsory parameter is missing an error will be raised at
+start-up.
+
+Finally, SWIFT outputs two YAML files at the start of a run. The first
+one ``used_parameters.yml`` contains all the parameters that were used
+for this run, **including all the optional parameters with their
+default values**. This file can be used to start an exact copy of the
+run. The second file, ``unused_parameters.yml`` contains all the
+values that were not read from the parameter file. This can be used to
+simplify the parameter file or check that nothing important was
+ignored (for instance because the code is not configured to use some
+options).
+
+The rest of this page describes all the SWIFT parameters, split by
+section. A list of all the possible parameters is kept in the file
+``examples/parameter_examples.yml``.
+
+Internal Unit System
+--------------------
+
+This section describes the units used internally by the code. This is
+the system of units in which all the equations are solved. All
+physical constants are converted to this system and if the ICs use a
+different system (see :ref:`ICs_units_label`) the particle quantities
+will be converted when read in.
+
+The system of units is described using the value of the 5 basic units
+of any system with respect to the CGS system. Instead of using a unit
+of time we use a unit of velocity as this is more intuitive. Users
+hence need to provide:
+
+* a unit of length: ``UnitLength_in_cgs``,
+* a unit of mass: ``UnitMass_in_cgs``,
+* a unit of velocity ``UnitVelocity_in_cgs``,
+* a unit of electric current ``UnitCurrent_in_cgs``,
+* a unit of temperature ``UnitTemp_in_cgs``.
+
+All these need to be expressed with respect to their cgs counter-part
+(i.e. :math:`cm`, :math:`g`, :math:`cm/s`, :math:`A` and :math:`K`). Recall
+that there are no h-factors in any of SWIFT's quantities; we, for instance,
+use :math:`cm` and not :math:`cm/h`.
+
+For instance to use the commonly adopted system of 10^10 Msun as a
+unit for mass, mega-parsec as a unit of length and km/s as a unit of
+speed, we would use:
+
+.. code:: YAML
+
+   # Common unit system for cosmo sims
+   InternalUnitSystem:
+     UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+     UnitLength_in_cgs:   3.08567758e24 # 1 Mpc in centimeters
+     UnitVelocity_in_cgs: 1e5           # 1 km/s in centimeters per second
+     UnitCurrent_in_cgs:  1             # 1 Ampere
+     UnitTemp_in_cgs:     1             # 1 Kelvin   
+	  
+Note that there are currently no variables in any of the SWIFT physics
+schemes that make use of the unit of electric current. There is also
+no incentive to use anything else than Kelvin but that makes the whole
+system consistent with any possible unit system.
+
+If one is interested in using the more humourous `FFF unit
+system <https://en.wikipedia.org/wiki/FFF_system>`_ one would use
+
+.. code:: YAML
+
+   # FFF unit system
+   InternalUnitSystem:
+     UnitMass_in_cgs:     40823.3133  # 1 Firkin (fir) in grams
+     UnitLength_in_cgs:   20116.8     # 1 Furlong (fur) in cm
+     UnitVelocity_in_cgs: 0.01663095  # 1 Furlong (fur) per Fortnight (ftn) in cm/s
+     UnitCurrent_in_cgs:  1           # 1 Ampere
+     UnitTemp_in_cgs:     1           # 1 Kelvin   
+
+The value of the physical constants in this system is left as an
+exercise for the reader [#f1]_.
+
+Cosmology
+---------
+
+When running a cosmological simulation, this section set the values of the
+cosmological model. The epanded :math:`\Lambda\rm{CDM}` parameters governing the
+background evolution of the Univese need to be specified here. These are:
+
+* The reduced Hubble constant: :math:`h`: ``h``,
+* The matter density parameter :math:`\Omega_m`: ``Omega_m``,
+* The cosmological constant density parameter :math:`\Omega_\Lambda`: ``Omega_lambda``,
+* The baryon density parameter :math:`\Omega_b`: ``Omega_b``,
+* The radiation density parameter :math:`\Omega_r`: ``Omega_r``.
+
+The last parameter can be omitted and will default to :math:`\Omega_r = 0`.
+
+This section als specifies the start and end of the simulation expressed in
+terms of scale-factors. The two parameters are:
+
+* Initial scale-factor: ``a_begin``,
+* Final scale-factor: ``a_end``.
+
+Two additional optional parameters can be used to change the equation of
+state of dark energy :math:`w(a)`. We use the evolution law :math:`w(a) =
+w_0 + w_a (1 - a)`. The two parameters in the YAML file are:
+
+* The :math:`z=0` dark energy equation of state parameter :math:`w_0`: ``w_0``
+* The dark energy equation of state evolutio parameter :math:`w_a`: ``w_a``
+
+If unspecified these parameters default to the default
+:math:`\Lambda\rm{CDM}` values of :math:`w_0 = -1` and :math:`w_a = 0`.
+
+For a Planck+13 cosmological model (ignoring radiation density as is
+commonly done) and running from :math:`z=127` to :math:`z=0`, one would hence
+use the following parameters:
+
+.. code:: YAML
+
+   Cosmology:
+     a_begin:        0.0078125     # z = 127
+     a_end:          1.0           # z = 0
+     h:              0.6777        
+     Omega_m:        0.307         
+     Omega_lambda:   0.693         
+     Omega_b:        0.0455        
+     Omega_r:        0.            # (Optional)
+     w_0:            -1.0          # (Optional)
+     w_a:            0.            # (Optional)
+
+When running a non-cosmological simulation (i.e. without the ``-c`` runtime
+flag) this section of the YAML file is entirely ignored.
+     
+Gravity
+-------
+
+SPH
+---
+
+Time Integration
+----------------
+
+Physical Constants
+------------------
+
+For some idealised test it can be useful to overwrite the value of
+some physical constants. In particular the value of the gravitational
+constant. SWIFT offers an optional parameter to overwrite the value of
+that constant.
+
+.. code:: YAML
+
+   PhysicalConstants:
+     G:   1
+
+Note that this set :math:`G` to the specified value in the internal system
+of units. Setting a value of `1` when using the system of units (10^10 Msun,
+Mpc, km/s) will mean that :math:`G_N=1` in these units [#f2]_ instead of the
+normal value :math:`G_N=43.00927`.
+
+This option is only used for specific tests and debugging. This entire
+section of the YAML file can typically be left out.
+
+Snapshots
+---------
+
+Some additional specific options for the snapshot outputs are described in the
+following pages:
+
+.. toctree::
+   :maxdepth: 1
+
+   output_selection
+
+Statistics
+----------
+
+Restarts
+--------
+
+Scheduler
+---------
+
+Domain Decomposition
+--------------------
+
+.. [#f1] The thorough reader (or overly keen SWIFT tester) would find  that the speed of light is :math:`c=1.8026\times10^{12}\,\rm{fur}\,\rm{ftn}^{-1}`, Newton's contant becomes :math:`G_N=4.896735\times10^{-4}~\rm{fur}^3\,\rm{fir}^{-1}\,\rm{ftn}^{-2}` and Planck's constant turns into :math:`h=4.851453\times 10^{-34}~\rm{fur}^2\,\rm{fir}\,\rm{ftn}^{-1}`.
+
+
+.. [#f2] which would translate into a constant :math:`G_N=1.5517771\times10^{-9}~cm^{3}\,g^{-1}\,s^{-2}` if expressed in the CGS system.
diff --git a/doc/RTD/source/GettingStarted/parameter_file.rst b/doc/RTD/source/ParameterFiles/output_selection.rst
similarity index 69%
rename from doc/RTD/source/GettingStarted/parameter_file.rst
rename to doc/RTD/source/ParameterFiles/output_selection.rst
index 550040ed25ec307633d6fade81eced58ed65a254..ca905c5e613082cdf9de9db718bdd16c4a3c8951 100644
--- a/doc/RTD/source/GettingStarted/parameter_file.rst
+++ b/doc/RTD/source/ParameterFiles/output_selection.rst
@@ -1,23 +1,17 @@
 .. Parameter File
    Loic Hausammann, 1 june 2018
 
-.. _Parameter_File_label:
-
-Parameter File
-==============
-
-To run SWIFT, you will need to provide a ``yaml`` parameter file.  An example is
-given in ``examples/parameter_file.yml`` which should contain all possible
-parameters.  Each section in this file corresponds to a different option in
-SWIFT and are not always required depending on the configuration options and
-the run time parameters.
+.. _Output_list_label:
 
 Output List
 ~~~~~~~~~~~
 
-In the sections ``Snapshots`` and ``Statistics``, you can specify the options ``output_list_on`` and ``output_list``  which receive an int and a filename.
-The ``output_list_on`` enable or not the output list and ``output_list`` is the filename containing the output times.
-With the file header, you can choose between writing redshifts, scale factors or times.
+In the sections ``Snapshots`` and ``Statistics``, you can specify the
+options ``output_list_on`` and ``output_list`` which receive an int
+and a filename.  The ``output_list_on`` enable or not the output list
+and ``output_list`` is the filename containing the output times.  With
+the file header, you can choose between writing redshifts, scale
+factors or times.
 
 Example of file containing with times (in internal units)::
 
@@ -42,6 +36,8 @@ Example of file with redshift::
   10
   5
 
+.. _Output_selection_label:
+  
 Output Selection
 ~~~~~~~~~~~~~~~~
 
diff --git a/doc/RTD/source/index.rst b/doc/RTD/source/index.rst
index 21eb7253c62c09595dd547540d9923c1be866e6a..343713f5f8c466013ea7f4ffaca0f010827dc404 100644
--- a/doc/RTD/source/index.rst
+++ b/doc/RTD/source/index.rst
@@ -15,6 +15,8 @@ difference is the parameter file that will need to be adapted for SWIFT.
    :maxdepth: 2
 
    GettingStarted/index
+   CommandLineOptions/index
+   ParameterFiles/index
    InitialConditions/index
    HydroSchemes/index
    Cooling/index
diff --git a/examples/DwarfGalaxy/README b/examples/DwarfGalaxy/README
new file mode 100644
index 0000000000000000000000000000000000000000..7a9167694a24c088997316180233b28b9126f298
--- /dev/null
+++ b/examples/DwarfGalaxy/README
@@ -0,0 +1,7 @@
+This example is a galaxy extracted from the example "ZoomIn". It allows
+to test SWIFT on a smaller problem. See the README in "ZoomIn" for more
+information.
+
+
+MD5 check-sum of the ICS: 
+ae2af84d88f30011b6a8af3f37d140cf  dwarf_galaxy.hdf5
\ No newline at end of file
diff --git a/examples/DwarfGalaxy/dwarf_galaxy.yml b/examples/DwarfGalaxy/dwarf_galaxy.yml
new file mode 100644
index 0000000000000000000000000000000000000000..655d41a92827ee56c08e977fe895d28ea8541d9e
--- /dev/null
+++ b/examples/DwarfGalaxy/dwarf_galaxy.yml
@@ -0,0 +1,71 @@
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e21 # kpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+# Structure finding options
+StructureFinding:
+  config_file_name:     stf_input.cfg # Name of the STF config file.
+  basename:             ./stf         # Common part of the name of output files.
+  output_time_format:   0             # Specifies the frequency format of structure finding. 0 for simulation steps (delta_step) and 1 for simulation time intervals (delta_time).
+  scale_factor_first:   0.92          # Scale-factor of the first snaphot (cosmological run)
+  time_first:           0.01        # Time of the first structure finding output (in internal units).
+  delta_step:           1000          # Time difference between consecutive structure finding outputs (in internal units) in simulation steps.
+  delta_time:           1.10          # Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
+
+# Cosmological parameters
+Cosmology:
+  h:              0.673        # Reduced Hubble constant
+  a_begin:        0.9873046739     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.315         # Matter density parameter
+  Omega_lambda:   0.685         # Dark-energy density parameter
+  Omega_b:        0.0486        # Baryon density parameter
+  
+Scheduler:
+  max_top_level_cells:    8
+  cell_split_size:           400       # (Optional) Maximal number of particles per cell (this is the default value).
+  
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            dwarf_galaxy # Common part of the name of output files
+  time_first:          0.  # Time of the first output (non-cosmological run) (in internal units)
+  delta_time:          0.02  # Time difference between consecutive outputs (in internal units)
+  compression:         1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.987345 # Scale-factor of the first stat dump (cosmological run)
+  time_first:          0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
+  delta_time:          1.05 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                    0.025    # Constant dimensionless multiplier for time integration.
+  theta:                  0.7      # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.05 # Comoving softening length (in internal units).
+  max_physical_softening: 0.01    # Physical softening length (in internal units).
+  mesh_side_length:       16
+
+# Parameters for the hydrodynamics scheme
+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.
+  minimal_temperature:   100      # (internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./dwarf_galaxy.hdf5     # The file to read
+  cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
+  cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
+
diff --git a/examples/DwarfGalaxy/getIC.sh b/examples/DwarfGalaxy/getIC.sh
new file mode 100644
index 0000000000000000000000000000000000000000..92f4cd3939845d57a61683e95135163b8639371f
--- /dev/null
+++ b/examples/DwarfGalaxy/getIC.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget https://obswww.unige.ch/~lhausamm/swift/IC/DwarfGalaxy/dwarf_galaxy.hdf5
diff --git a/examples/DwarfGalaxy/run.sh b/examples/DwarfGalaxy/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..4cc87880215a37c9eac59b19e584f4887cba2c38
--- /dev/null
+++ b/examples/DwarfGalaxy/run.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e dwarf_galaxy.hdf5 ]
+then
+    echo "Fetching initial conditions for the dwarf galaxy example..."
+    ./getIC.sh
+fi
+
+../swift -b -G -s -S -t 8 dwarf_galaxy.yml 2>&1 | tee output.log
+
diff --git a/examples/ZoomIn/README b/examples/ZoomIn/README
new file mode 100644
index 0000000000000000000000000000000000000000..cffc275f2ae1046156d392f8725a7b542c80471a
--- /dev/null
+++ b/examples/ZoomIn/README
@@ -0,0 +1,16 @@
+Initial conditions for a zoom in cosmological simulation of dwarf
+galaxies. These have been generated by MUSIC and ran up to z=0 with
+GEAR (see Revaz and Jablonka 2018 for more details on the simulation).
+
+The cosmology is taken from Planck 2015.
+
+The initial conditions have been cleaned to contain only the required
+fields. The ICs have been created for Gadget and the positions and box
+size are hence expressed in h-full units (e.g. box size of 32 / h Mpc).
+Similarly, the peculiar velocitites contain an extra sqrt(a) factor. 
+
+We will use SWIFT to cancel the h- and a-factors from the ICs. Gas
+particles will be generated at startup.
+
+MD5 check-sum of the ICS: 
+9aafe154438478ed435e88664c1c5dba zoom_in.hdf5
diff --git a/examples/ZoomIn/getIC.sh b/examples/ZoomIn/getIC.sh
new file mode 100644
index 0000000000000000000000000000000000000000..6cdfaec981af515249578faa72798c53448e7ecb
--- /dev/null
+++ b/examples/ZoomIn/getIC.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget https://obswww.unige.ch/~lhausamm/swift/IC/ZoomIn/zoom_in.hdf5
diff --git a/examples/ZoomIn/run.sh b/examples/ZoomIn/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..99eda1cfc779c1958d19d0c7ae234b6c211f8127
--- /dev/null
+++ b/examples/ZoomIn/run.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e zoom_in.hdf5 ]
+then
+    echo "Fetching initial conditions for the zoom in example..."
+    ./getIC.sh
+fi
+
+../swift -b -c -G -s -S -t 8 zoom_in.yml 2>&1 | tee output.log
+
diff --git a/examples/ZoomIn/zoom_in.yml b/examples/ZoomIn/zoom_in.yml
new file mode 100644
index 0000000000000000000000000000000000000000..5f0fa4958b63150d94ff06df042a03af6f221038
--- /dev/null
+++ b/examples/ZoomIn/zoom_in.yml
@@ -0,0 +1,61 @@
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e21 # kpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+# Cosmological parameters
+Cosmology:
+  h:              0.673        # Reduced Hubble constant
+  a_begin:        0.9873046739     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.315         # Matter density parameter
+  Omega_lambda:   0.685         # Dark-energy density parameter
+  Omega_b:        0.0486        # Baryon density parameter
+  
+Scheduler:
+  max_top_level_cells:    8
+  
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1e-2  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            zoom_in # Common part of the name of output files
+  scale_factor_first:  0.987345  # Scale-factor of the first snaphot (cosmological run)
+  time_first:          0.01  # Time of the first output (non-cosmological run) (in internal units)
+  delta_time:          1.01  # Time difference between consecutive outputs (in internal units)
+  compression:         1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.987345 # Scale-factor of the first stat dump (cosmological run)
+  time_first:          0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
+  delta_time:          1.05 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                    0.025    # Constant dimensionless multiplier for time integration.
+  theta:                  0.7      # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.05 # Comoving softening length (in internal units).
+  max_physical_softening: 0.01    # Physical softening length (in internal units).
+  mesh_side_length:       16
+
+# Parameters for the hydrodynamics scheme
+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.
+  minimal_temperature:   100      # (internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./zoom_in.hdf5     # The file to read
+  cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
+  cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
+
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index ab422fd8f2f29a6b9fd200d4cfd6ea26ca120497..f86abe054e596be123328497bbb340cd30abf3e8 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -34,6 +34,10 @@ SPH:
   minimal_temperature:   0        # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
   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
 
 # Parameters for the self-gravity scheme
 Gravity:
@@ -133,9 +137,9 @@ DomainDecomposition:
                             # new decomposition, or number of steps (>1) between decompositions
   minfrac:          0.9     # (Optional) Fractional of all particles that should be updated in previous step when
                             # using CPU time trigger
-  usemetis          0       # Use serial METIS when ParMETIS is also available.
-  adaptive          1       # Use adaptive repartition when ParMETIS is available, otherwise simple refinement.
-  itr               100     # When adaptive defines the ratio of inter node communication time to data redistribution time, in the range 0.00001 to 10000000.0.
+  usemetis:         0       # Use serial METIS when ParMETIS is also available.
+  adaptive:         1       # Use adaptive repartition when ParMETIS is available, otherwise simple refinement.
+  itr:              100     # When adaptive defines the ratio of inter node communication time to data redistribution time, in the range 0.00001 to 10000000.0.
                             # Lower values give less data movement during redistributions, at the cost of global balance which may require more communication.
 
 # Parameters related to the equation of state ------------------------------------------
diff --git a/examples/plot_task_level.py b/examples/plot_task_level.py
new file mode 100644
index 0000000000000000000000000000000000000000..2d51f4829692bc71826ec6c4b93c025f7106828f
--- /dev/null
+++ b/examples/plot_task_level.py
@@ -0,0 +1,43 @@
+#!/usr/bin/python
+"""
+Usage:
+  ./plot_task_level.py task_level.txt
+
+Description:
+  Plot the number of tasks for each depth level and each type of task.
+"""
+
+
+import pandas as pd
+import matplotlib.pyplot as plt
+import sys
+
+# get filename
+filename = sys.argv[-1]
+
+# Column names
+names = ["type", "subtype", "depth", "count"]
+
+# read file
+data = pd.read_csv(filename, sep=' ', comment="#", names=names)
+
+# generate color map
+cmap = plt.get_cmap("hsv")
+N = data["depth"].max() + 5
+
+# plot data
+for i in range(data["depth"].max()):
+    ind = data["depth"] == i
+    label = "depth = %i" % i
+    c = cmap(i / N)
+    plt.plot(data["type"][ind] + "_" + data["subtype"][ind],
+             data["count"][ind], ".", label=label,
+             color=c)
+
+# modify figure parameters and show it
+plt.gca().set_yscale("log")
+plt.xticks(rotation=45)
+plt.ylabel("Number of Tasks")
+plt.gcf().subplots_adjust(bottom=0.15)
+plt.legend()
+plt.show()
diff --git a/examples/plot_tasks.py b/examples/plot_tasks.py
index ec184ca1feee5879366d31ef8a984879bd0cddaf..42d362bfd3cd11a68c6c4349af4e21ef100e1578 100755
--- a/examples/plot_tasks.py
+++ b/examples/plot_tasks.py
@@ -157,11 +157,11 @@ for task in SUBTYPES:
 
 #  For fiddling with colours...
 if args.verbose:
-    print "#Selected colours:"
+    print("#Selected colours:")
     for task in sorted(TASKCOLOURS.keys()):
-        print "# " + task + ": " + TASKCOLOURS[task]
+        print("# " + task + ": " + TASKCOLOURS[task])
     for task in sorted(SUBCOLOURS.keys()):
-        print "# " + task + ": " + SUBCOLOURS[task]
+        print("# " + task + ": " + SUBCOLOURS[task])
 
 #  Read input.
 data = pl.loadtxt( infile )
@@ -169,11 +169,11 @@ data = pl.loadtxt( infile )
 #  Do we have an MPI file?
 full_step = data[0,:]
 if full_step.size == 13:
-    print "# MPI mode"
+    print("# MPI mode")
     mpimode = True
     if ranks == None:
         ranks = range(int(max(data[:,0])) + 1)
-    print "# Number of ranks:", len(ranks)
+    print("# Number of ranks:", len(ranks))
     rankcol = 0
     threadscol = 1
     taskcol = 2
@@ -181,7 +181,7 @@ if full_step.size == 13:
     ticcol = 5
     toccol = 6
 else:
-    print "# non MPI mode"
+    print("# non MPI mode")
     ranks = [0]
     mpimode = False
     rankcol = -1
@@ -194,10 +194,10 @@ else:
 #  Get CPU_CLOCK to convert ticks into milliseconds.
 CPU_CLOCK = float(full_step[-1]) / 1000.0
 if args.verbose:
-    print "# CPU frequency:", CPU_CLOCK * 1000.0
+    print("# CPU frequency:", CPU_CLOCK * 1000.0)
 
 nthread = int(max(data[:,threadscol])) + 1
-print "# Number of threads:", nthread
+print("# Number of threads:", nthread)
 
 # Avoid start and end times of zero.
 sdata = data[data[:,ticcol] != 0]
@@ -224,24 +224,24 @@ if delta_t == 0:
         dt = toc_step - tic_step
         if dt > delta_t:
             delta_t = dt
-    print "# Data range: ", delta_t / CPU_CLOCK, "ms"
+    print("# Data range: ", delta_t / CPU_CLOCK, "ms")
 
 # Once more doing the real gather and plots this time.
 for rank in ranks:
-    print "# Processing rank: ", rank
+    print("# Processing rank: ", rank)
     if mpimode:
         data = sdata[sdata[:,rankcol] == rank]
         full_step = data[0,:]
     tic_step = int(full_step[ticcol])
     toc_step = int(full_step[toccol])
-    print "# Min tic = ", tic_step
+    print("# Min tic = ", tic_step)
     data = data[1:,:]
     typesseen = []
     nethread = 0
 
     #  Dummy image for ranks that have no tasks.
     if data.size == 0:
-        print "# Rank ", rank, " has no tasks"
+        print("# Rank ", rank, " has no tasks")
         fig = pl.figure()
         ax = fig.add_subplot(1,1,1)
         ax.set_xlim(-delta_t * 0.01 / CPU_CLOCK, delta_t * 1.01 / CPU_CLOCK)
@@ -367,6 +367,6 @@ for rank in ranks:
     else:
         outpng = outbase + ".png"
     pl.savefig(outpng)
-    print "Graphics done, output written to", outpng
+    print("Graphics done, output written to", outpng)
 
 sys.exit(0)
diff --git a/src/cache.h b/src/cache.h
index 77f933cc54e7549ef55b7e7503b9797524413713..5dd8164b1dc80795a8593cc2af42c2c9e7e68885 100644
--- a/src/cache.h
+++ b/src/cache.h
@@ -834,6 +834,7 @@ static INLINE void cache_clean(struct cache *c) {
     free(c->balsara);
     free(c->soundspeed);
   }
+  c->count = 0;
 }
 
 #endif /* WITH_VECTORIZATION */
diff --git a/src/cell.c b/src/cell.c
index b17ca11f2b09650407a164fc2cc09428aadbe353..60e2e9393fe62f381b598a807fc7a631b7d8dbd5 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -337,6 +337,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
       temp->split = 0;
       temp->hydro.dx_max_part = 0.f;
       temp->hydro.dx_max_sort = 0.f;
+      temp->stars.dx_max_part = 0.f;
       temp->nodeID = c->nodeID;
       temp->parent = c;
       c->progeny[k] = temp;
@@ -410,6 +411,7 @@ int cell_pack_end_step(struct cell *restrict c,
   pcells[0].grav.ti_end_min = c->grav.ti_end_min;
   pcells[0].grav.ti_end_max = c->grav.ti_end_max;
   pcells[0].hydro.dx_max_part = c->hydro.dx_max_part;
+  pcells[0].stars.dx_max_part = c->stars.dx_max_part;
 
   /* Fill in the progeny, depth-first recursion. */
   int count = 1;
@@ -446,6 +448,7 @@ int cell_unpack_end_step(struct cell *restrict c,
   c->grav.ti_end_min = pcells[0].grav.ti_end_min;
   c->grav.ti_end_max = pcells[0].grav.ti_end_max;
   c->hydro.dx_max_part = pcells[0].hydro.dx_max_part;
+  c->stars.dx_max_part = pcells[0].stars.dx_max_part;
 
   /* Fill in the progeny, depth-first recursion. */
   int count = 1;
@@ -1156,8 +1159,11 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
 void cell_sanitize(struct cell *c, int treated) {
 
   const int count = c->hydro.count;
+  const int scount = c->stars.count;
   struct part *parts = c->hydro.parts;
+  struct spart *sparts = c->stars.parts;
   float h_max = 0.f;
+  float stars_h_max = 0.f;
 
   /* Treat cells will <1000 particles */
   if (count < 1000 && !treated) {
@@ -1170,6 +1176,10 @@ void cell_sanitize(struct cell *c, int treated) {
       if (parts[i].h == 0.f || parts[i].h > upper_h_max)
         parts[i].h = upper_h_max;
     }
+    for (int i = 0; i < scount; ++i) {
+      if (sparts[i].h == 0.f || sparts[i].h > upper_h_max)
+        sparts[i].h = upper_h_max;
+    }
   }
 
   /* Recurse and gather the new h_max values */
@@ -1183,16 +1193,20 @@ void cell_sanitize(struct cell *c, int treated) {
 
         /* And collect */
         h_max = max(h_max, c->progeny[k]->hydro.h_max);
+        stars_h_max = max(stars_h_max, c->progeny[k]->stars.h_max);
       }
     }
   } else {
 
     /* Get the new value of h_max */
     for (int i = 0; i < count; ++i) h_max = max(h_max, parts[i].h);
+    for (int i = 0; i < scount; ++i)
+      stars_h_max = max(stars_h_max, sparts[i].h);
   }
 
   /* Record the change */
   c->hydro.h_max = h_max;
+  c->stars.h_max = stars_h_max;
 }
 
 /**
@@ -1624,6 +1638,13 @@ void cell_activate_drift_gpart(struct cell *c, struct scheduler *s) {
   }
 }
 
+/**
+ * @brief Activate the #spart drifts on the given cell.
+ */
+void cell_activate_drift_spart(struct cell *c, struct scheduler *s) {
+  cell_activate_drift_gpart(c, s);
+}
+
 /**
  * @brief Activate the sorts up a cell hierarchy.
  */
@@ -1695,6 +1716,7 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
   /* Store the current dx_max and h_max values. */
   ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
   ci->hydro.h_max_old = ci->hydro.h_max;
+
   if (cj != NULL) {
     cj->hydro.dx_max_part_old = cj->hydro.dx_max_part;
     cj->hydro.h_max_old = cj->hydro.h_max;
@@ -2045,7 +2067,353 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
  */
 void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
                                        struct scheduler *s) {
-  // LOIC: to implement
+  const struct engine *e = s->space->e;
+
+  /* Store the current dx_max and h_max values. */
+  ci->stars.dx_max_part_old = ci->stars.dx_max_part;
+  ci->stars.h_max_old = ci->stars.h_max;
+
+  if (cj != NULL) {
+    cj->stars.dx_max_part_old = cj->stars.dx_max_part;
+    cj->stars.h_max_old = cj->stars.h_max;
+  }
+
+  /* Self interaction? */
+  if (cj == NULL) {
+
+    /* Do anything? */
+    if (ci->stars.count == 0 || !cell_is_active_stars(ci, e)) return;
+
+    /* Recurse? */
+    if (cell_can_recurse_in_self_stars_task(ci)) {
+
+      /* Loop over all progenies and pairs of progenies */
+      for (int j = 0; j < 8; j++) {
+        if (ci->progeny[j] != NULL) {
+          cell_activate_subcell_stars_tasks(ci->progeny[j], NULL, s);
+          for (int k = j + 1; k < 8; k++)
+            if (ci->progeny[k] != NULL)
+              cell_activate_subcell_stars_tasks(ci->progeny[j], ci->progeny[k],
+                                                s);
+        }
+      }
+    } else {
+
+      /* We have reached the bottom of the tree: activate drift */
+      cell_activate_drift_spart(ci, s);
+      cell_activate_drift_part(ci, s);
+    }
+  }
+
+  /* Otherwise, pair interation */
+  else {
+
+    /* Should we even bother? */
+    if (!cell_is_active_stars(ci, e) && !cell_is_active_stars(cj, e)) return;
+    if (ci->stars.count == 0 || cj->stars.count == 0) return;
+
+    /* Get the orientation of the pair. */
+    double shift[3];
+    int sid = space_getsid(s->space, &ci, &cj, shift);
+
+    /* recurse? */
+    if (cell_can_recurse_in_pair_stars_task(ci) &&
+        cell_can_recurse_in_pair_stars_task(cj)) {
+
+      /* Different types of flags. */
+      switch (sid) {
+
+        /* Regular sub-cell interactions of a single cell. */
+        case 0: /* (  1 ,  1 ,  1 ) */
+          if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[0],
+                                              s);
+          break;
+
+        case 1: /* (  1 ,  1 ,  0 ) */
+          if (ci->progeny[6] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[0],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[1],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[0],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[1],
+                                              s);
+          break;
+
+        case 2: /* (  1 ,  1 , -1 ) */
+          if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[1],
+                                              s);
+          break;
+
+        case 3: /* (  1 ,  0 ,  1 ) */
+          if (ci->progeny[5] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[0],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[2],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[0],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[2],
+                                              s);
+          break;
+
+        case 4: /* (  1 ,  0 ,  0 ) */
+          if (ci->progeny[4] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[4], cj->progeny[0],
+                                              s);
+          if (ci->progeny[4] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[4], cj->progeny[1],
+                                              s);
+          if (ci->progeny[4] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[4], cj->progeny[2],
+                                              s);
+          if (ci->progeny[4] != NULL && cj->progeny[3] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[4], cj->progeny[3],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[0],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[1],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[2],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[3] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[3],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[0],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[1],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[2],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[3] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[3],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[0],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[1],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[2],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[3] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[3],
+                                              s);
+          break;
+
+        case 5: /* (  1 ,  0 , -1 ) */
+          if (ci->progeny[4] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[4], cj->progeny[1],
+                                              s);
+          if (ci->progeny[4] != NULL && cj->progeny[3] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[4], cj->progeny[3],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[1],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[3] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[3],
+                                              s);
+          break;
+
+        case 6: /* (  1 , -1 ,  1 ) */
+          if (ci->progeny[5] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[2],
+                                              s);
+          break;
+
+        case 7: /* (  1 , -1 ,  0 ) */
+          if (ci->progeny[4] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[4], cj->progeny[2],
+                                              s);
+          if (ci->progeny[4] != NULL && cj->progeny[3] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[4], cj->progeny[3],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[2],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[3] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[3],
+                                              s);
+          break;
+
+        case 8: /* (  1 , -1 , -1 ) */
+          if (ci->progeny[4] != NULL && cj->progeny[3] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[4], cj->progeny[3],
+                                              s);
+          break;
+
+        case 9: /* (  0 ,  1 ,  1 ) */
+          if (ci->progeny[3] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[3], cj->progeny[0],
+                                              s);
+          if (ci->progeny[3] != NULL && cj->progeny[4] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[3], cj->progeny[4],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[0],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[4] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[4],
+                                              s);
+          break;
+
+        case 10: /* (  0 ,  1 ,  0 ) */
+          if (ci->progeny[2] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[2], cj->progeny[0],
+                                              s);
+          if (ci->progeny[2] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[2], cj->progeny[1],
+                                              s);
+          if (ci->progeny[2] != NULL && cj->progeny[4] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[2], cj->progeny[4],
+                                              s);
+          if (ci->progeny[2] != NULL && cj->progeny[5] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[2], cj->progeny[5],
+                                              s);
+          if (ci->progeny[3] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[3], cj->progeny[0],
+                                              s);
+          if (ci->progeny[3] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[3], cj->progeny[1],
+                                              s);
+          if (ci->progeny[3] != NULL && cj->progeny[4] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[3], cj->progeny[4],
+                                              s);
+          if (ci->progeny[3] != NULL && cj->progeny[5] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[3], cj->progeny[5],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[0],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[1],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[4] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[4],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[5] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[5],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[0],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[1],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[4] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[4],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[5] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[5],
+                                              s);
+          break;
+
+        case 11: /* (  0 ,  1 , -1 ) */
+          if (ci->progeny[2] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[2], cj->progeny[1],
+                                              s);
+          if (ci->progeny[2] != NULL && cj->progeny[5] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[2], cj->progeny[5],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[1],
+                                              s);
+          if (ci->progeny[6] != NULL && cj->progeny[5] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[6], cj->progeny[5],
+                                              s);
+          break;
+
+        case 12: /* (  0 ,  0 ,  1 ) */
+          if (ci->progeny[1] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[1], cj->progeny[0],
+                                              s);
+          if (ci->progeny[1] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[1], cj->progeny[2],
+                                              s);
+          if (ci->progeny[1] != NULL && cj->progeny[4] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[1], cj->progeny[4],
+                                              s);
+          if (ci->progeny[1] != NULL && cj->progeny[6] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[1], cj->progeny[6],
+                                              s);
+          if (ci->progeny[3] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[3], cj->progeny[0],
+                                              s);
+          if (ci->progeny[3] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[3], cj->progeny[2],
+                                              s);
+          if (ci->progeny[3] != NULL && cj->progeny[4] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[3], cj->progeny[4],
+                                              s);
+          if (ci->progeny[3] != NULL && cj->progeny[6] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[3], cj->progeny[6],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[0],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[2],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[4] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[4],
+                                              s);
+          if (ci->progeny[5] != NULL && cj->progeny[6] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[5], cj->progeny[6],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[0],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[2] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[2],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[4] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[4],
+                                              s);
+          if (ci->progeny[7] != NULL && cj->progeny[6] != NULL)
+            cell_activate_subcell_stars_tasks(ci->progeny[7], cj->progeny[6],
+                                              s);
+          break;
+      }
+
+    }
+
+    /* Otherwise, activate the sorts and drifts. */
+    else if (cell_is_active_stars(ci, e) || cell_is_active_stars(cj, e)) {
+
+      /* We are going to interact this pair, so store some values. */
+      atomic_or(&ci->hydro.requires_sorts, 1 << sid);
+      atomic_or(&cj->hydro.requires_sorts, 1 << sid);
+      ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
+      cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
+
+      /* Activate the drifts if the cells are local. */
+      if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
+      if (ci->nodeID == engine_rank) cell_activate_drift_spart(ci, s);
+      if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
+      if (cj->nodeID == engine_rank) cell_activate_drift_spart(cj, s);
+
+      /* Do we need to sort the cells? */
+      cell_activate_sorts(ci, sid, s);
+      cell_activate_sorts(cj, sid, s);
+    }
+  } /* Otherwise, pair interation */
 }
 
 /**
@@ -2999,11 +3367,15 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
  */
 void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
 
+  const float stars_h_max = e->stars_properties->h_max;
   const integertime_t ti_old_gpart = c->grav.ti_old_part;
   const integertime_t ti_current = e->ti_current;
   struct gpart *const gparts = c->grav.parts;
   struct spart *const sparts = c->stars.parts;
 
+  float dx_max = 0.f, dx2_max = 0.f;
+  float cell_h_max = 0.f;
+
   /* Drift irrespective of cell flags? */
   force |= c->grav.do_drift;
 
@@ -3025,9 +3397,17 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
 
         /* Recurse */
         cell_drift_gpart(cp, e, force);
+
+        /* Update */
+        dx_max = max(dx_max, cp->stars.dx_max_part);
+        cell_h_max = max(cell_h_max, cp->stars.h_max);
       }
     }
 
+    /* Store the values */
+    c->stars.h_max = cell_h_max;
+    c->stars.dx_max_part = dx_max;
+
     /* Update the time of the last drift */
     c->grav.ti_old_part = ti_current;
 
@@ -3111,8 +3491,26 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
       }
 #endif
 
-      /* Note: no need to compute dx_max as all spart have a gpart */
-    }
+      /* Limit h to within the allowed range */
+      sp->h = min(sp->h, stars_h_max);
+
+      /* Compute (square of) motion since last cell construction */
+      const float dx2 = sp->x_diff[0] * sp->x_diff[0] +
+                        sp->x_diff[1] * sp->x_diff[1] +
+                        sp->x_diff[2] * sp->x_diff[2];
+      dx2_max = max(dx2_max, dx2);
+
+      /* Maximal smoothing length */
+      cell_h_max = max(cell_h_max, sp->h);
+
+    } /* Note: no need to compute dx_max as all spart have a gpart */
+
+    /* Now, get the maximal particle motion from its square */
+    dx_max = sqrtf(dx2_max);
+
+    /* Store the values */
+    c->stars.h_max = cell_h_max;
+    c->stars.dx_max_part = dx_max;
 
     /* Update the time of the last drift */
     c->grav.ti_old_part = ti_current;
diff --git a/src/cell.h b/src/cell.h
index e512da5473d70e55c7937cf42eca65d7e82a06f5..a01f2eaaaf868107ac3e0fddd1540962158913c3 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -144,6 +144,9 @@ struct pcell {
     /*! Number of #spart in this cell. */
     int count;
 
+    /*! Maximal smoothing length. */
+    double h_max;
+
   } stars;
 
   /*! Relative indices of the cell's progeny. */
@@ -188,6 +191,10 @@ struct pcell_step {
 
   /*! Stars variables */
   struct {
+
+    /*! Maximal distance any #part has travelled since last rebuild */
+    float dx_max_part;
+
   } stars;
 };
 
@@ -435,6 +442,18 @@ struct cell {
     /*! Nr of #spart in this cell. */
     int count;
 
+    /*! Max smoothing length in this cell. */
+    double h_max;
+
+    /*! Values of h_max before the drifts, used for sub-cell tasks. */
+    float h_max_old;
+
+    /*! Maximum part movement in this cell since last construction. */
+    float dx_max_part;
+
+    /*! Values of dx_max before the drifts, used for sub-cell tasks. */
+    float dx_max_part_old;
+
     /*! Dependency implicit task for the star ghost  (in->ghost->out)*/
     struct task *ghost_in;
 
@@ -684,8 +703,12 @@ cell_can_recurse_in_self_hydro_task(const struct cell *c) {
 __attribute__((always_inline)) INLINE static int
 cell_can_recurse_in_pair_stars_task(const struct cell *c) {
 
-  // LOIC: To implement
-  return 0;
+  /* Is the cell split ? */
+  /* If so, is the cut-off radius plus the max distance the parts have moved */
+  /* smaller than the sub-cell sizes ? */
+  /* Note: We use the _old values as these might have been updated by a drift */
+  return c->split && ((kernel_gamma * c->stars.h_max_old +
+                       c->stars.dx_max_part_old) < 0.5f * c->dmin);
 }
 
 /**
@@ -697,8 +720,8 @@ cell_can_recurse_in_pair_stars_task(const struct cell *c) {
 __attribute__((always_inline)) INLINE static int
 cell_can_recurse_in_self_stars_task(const struct cell *c) {
 
-  // LOIC: To implement
-  return 0;
+  /* Is the cell split and not smaller than the smoothing length? */
+  return c->split && (kernel_gamma * c->stars.h_max_old < 0.5f * c->dmin);
 }
 
 /**
@@ -746,8 +769,13 @@ __attribute__((always_inline)) INLINE static int cell_can_split_self_hydro_task(
 __attribute__((always_inline)) INLINE static int cell_can_split_pair_stars_task(
     const struct cell *c) {
 
-  // LOIC: To implement
-  return 0;
+  /* Is the cell split ? */
+  /* If so, is the cut-off radius with some leeway smaller than */
+  /* the sub-cell sizes ? */
+  /* Note that since tasks are create after a rebuild no need to take */
+  /* into account any part motion (i.e. dx_max == 0 here) */
+  return c->split &&
+         (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
 }
 
 /**
@@ -759,8 +787,13 @@ __attribute__((always_inline)) INLINE static int cell_can_split_pair_stars_task(
 __attribute__((always_inline)) INLINE static int cell_can_split_self_stars_task(
     const struct cell *c) {
 
-  // LOIC: To implement
-  return 0;
+  /* Is the cell split ? */
+  /* If so, is the cut-off radius with some leeway smaller than */
+  /* the sub-cell sizes ? */
+  /* Note: No need for more checks here as all the sub-pairs and sub-self */
+  /* tasks will be created. So no need to check for h_max */
+  return c->split &&
+         (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
 }
 
 /**
diff --git a/src/const.h b/src/const.h
index 821e671ac4ce61334b669f60a4785a7daf12d51d..e417b8ca3827ef87396706c56df36bb9bd3aed75 100644
--- a/src/const.h
+++ b/src/const.h
@@ -21,18 +21,11 @@
 #define SWIFT_CONST_H
 
 /* SPH Viscosity constants. */
-/* Cosmology defaults: a=0.8, b=3.0. Planetary defaults: a=1.5, b=4.0
+/* Cosmology default beta=3.0. Planetary default beta=4.0
+ * Alpha can be set in the parameter file.
  * Beta is defined as in e.g. Price (2010) Eqn (103) */
-#define const_viscosity_alpha 0.8f
 #define const_viscosity_beta 3.0f
 
-#define const_viscosity_alpha_min \
-  0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
-#define const_viscosity_alpha_max \
-  2.0f /* Values taken from (Price,2004), not used in legacy gadget mode */
-#define const_viscosity_length \
-  0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
-
 /* SPH Thermal conductivity constants. */
 #define const_conductivity_alpha \
   1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
diff --git a/src/drift.h b/src/drift.h
index ff0fea744012b7143afed2a05b286d4646cdd69a..351b15381dde9a95af908003b1c8df01b56610fa 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -137,6 +137,12 @@ __attribute__((always_inline)) INLINE static void drift_spart(
   sp->x[0] += sp->v[0] * dt_drift;
   sp->x[1] += sp->v[1] * dt_drift;
   sp->x[2] += sp->v[2] * dt_drift;
+
+  /* Compute offsets since last cell construction */
+  for (int k = 0; k < 3; k++) {
+    const float dx = sp->v[k] * dt_drift;
+    sp->x_diff[k] -= dx;
+  }
 }
 
 #endif /* SWIFT_DRIFT_H */
diff --git a/src/engine.c b/src/engine.c
index 2f09e9c9fbd3e5deb8d56c35f0793b2f37faa7a5..40365c4676c6ce8ad2d402b28fdea34ea606b597 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3379,14 +3379,11 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 }
 
 /**
- * @brief Duplicates the first hydro loop and construct all the
- * dependencies for the stars
+ * @brief Creates all the task dependencies for the stars
  *
- * This is done by looping over all the previously constructed tasks
- * and adding another task involving the same cells but this time
- * corresponding to the second hydro loop over neighbours.
- * With all the relevant tasks for a given cell available, we construct
- * all the dependencies for that cell.
+ * @param map_data The tasks
+ * @param num_elements number of tasks
+ * @param extra_data The #engine
  */
 void engine_link_stars_tasks_mapper(void *map_data, int num_elements,
                                     void *extra_data) {
@@ -3408,8 +3405,9 @@ void engine_link_stars_tasks_mapper(void *map_data, int num_elements,
 
       /* Now, build all the dependencies for the stars */
       engine_make_stars_loops_dependencies(sched, t, t->ci);
-      scheduler_addunlock(sched, t->ci->stars.ghost_out,
-                          t->ci->super->end_force);
+      if (t->ci == t->ci->super)
+        scheduler_addunlock(sched, t->ci->super->stars.ghost_out,
+                            t->ci->super->end_force);
     }
 
     /* Otherwise, pair interaction? */
@@ -5165,6 +5163,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
 #endif
 
   if (e->nodeID == 0) scheduler_write_dependencies(&e->sched, e->verbose);
+  if (e->nodeID == 0) scheduler_write_task_level(&e->sched);
 
   /* Run the 0th time-step */
   TIMER_TIC2;
@@ -5246,6 +5245,20 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
     }
   }
 
+  if (s->cells_top != NULL && s->nr_sparts > 0) {
+    for (int i = 0; i < s->nr_cells; i++) {
+      struct cell *c = &s->cells_top[i];
+      if (c->nodeID == engine_rank && c->stars.count > 0) {
+        float spart_h_max = c->stars.parts[0].h;
+        for (int k = 1; k < c->stars.count; k++) {
+          if (c->stars.parts[k].h > spart_h_max)
+            spart_h_max = c->stars.parts[k].h;
+        }
+        c->stars.h_max = max(spart_h_max, c->stars.h_max);
+      }
+    }
+  }
+
   clocks_gettime(&time2);
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -7525,18 +7538,11 @@ void engine_clean(struct engine *e) {
   }
   free(e->runners);
   free(e->snapshot_units);
-  if (e->output_list_snapshots) {
-    output_list_clean(e->output_list_snapshots);
-    free(e->output_list_snapshots);
-  }
-  if (e->output_list_stats) {
-    output_list_clean(e->output_list_stats);
-    free(e->output_list_stats);
-  }
-  if (e->output_list_stf) {
-    output_list_clean(e->output_list_stf);
-    free(e->output_list_stf);
-  }
+
+  output_list_clean(&e->output_list_snapshots);
+  output_list_clean(&e->output_list_stats);
+  output_list_clean(&e->output_list_stf);
+
   free(e->links);
   free(e->cell_loc);
   scheduler_clean(&e->sched);
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index b495af1d5b3ba1b91522bdeeac27fe77b262bd22..4252f2787aefcec058b8fa956eaa0351b8f41d57 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -458,12 +458,14 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @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 float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
   const float fac_mu = cosmo->a_factor_mu;
 
   /* Some smoothing length multiples. */
@@ -492,6 +494,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   p->force.balsara =
       normDiv_v / (normDiv_v + normRot_v + 0.0001f * fac_mu * fc * h_inv);
 
+  /* Set the AV property */
+  p->alpha = hydro_props->viscosity.alpha;
+
   /* Viscosity parameter decay time */
   /* const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
    */
@@ -500,8 +505,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   /* const float S = max(-normDiv_v, 0.f); */
 
   /* Compute the particle's viscosity parameter time derivative */
-  /* const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau + */
-  /*                         (const_viscosity_alpha_max - p->alpha) * S; */
+  /* const float alpha_dot = (hydro_props->viscosity.alpha_max) - p->alpha) /
+   * tau + */
+  /*                         (hydro_props->viscosity.alpha_max - p->alpha) * S;
+   */
 
   /* Update particle's viscosity paramter */
   /* p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase; */  // MATTHIEU
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index 51c448e7a6fa56216dc49197d8b9240ee5acf5f1..69919c202223fdecc197a87178e59767c02ee16e 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -183,13 +183,6 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
       h_grpsph, "Viscosity Model",
       "Morris & Monaghan (1997), Rosswog, Davies, Thielemann & "
       "Piran (2000) with additional Balsara (1995) switch");
-  io_write_attribute_f(h_grpsph, "Viscosity alpha_min",
-                       const_viscosity_alpha_min);
-  io_write_attribute_f(h_grpsph, "Viscosity alpha_max",
-                       const_viscosity_alpha_max);
-  io_write_attribute_f(h_grpsph, "Viscosity beta", 2.f);
-  io_write_attribute_f(h_grpsph, "Viscosity decay length",
-                       const_viscosity_length);
 
   /* Time integration properties */
   io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 36db326cf9624a77e3f2b4d06b5f67542f1edc9e..353fa6a021dfec59a7fd3fbe363e91b0495f4552 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -472,12 +472,14 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @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 float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
   const float fac_mu = cosmo->a_factor_mu;
 
@@ -502,8 +504,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float P_over_rho2 = pressure * rho_inv * rho_inv;
 
   /* Compute the Balsara switch */
+  /* Pre-multiply in the AV factor; hydro_props are not passed to the iact
+   * functions */
   const float balsara =
-      abs_div_v / (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
+      hydro_props->viscosity.alpha * abs_div_v /
+      (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
 
   /* Compute the "grad h" term */
   const float omega_inv =
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 05b74b6a88ca49969f14965dd55fee827627c4db..b7201520c0cd7b2c7039460890ca526ce25a2c02 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -497,8 +497,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -620,8 +619,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -654,8 +652,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 }
 
 #ifdef WITH_VECTORIZATION
-static const vector const_viscosity_alpha_fac =
-    FILL_VEC(-0.25f * const_viscosity_alpha);
 
 /**
  * @brief Force interaction computed using 1 vector
@@ -746,9 +742,9 @@ runner_iact_nonsym_1_vec_force(
 
   /* Now construct the full viscosity term */
   rho_ij.v = vec_mul(vec_set1(0.5f), vec_add(pirho.v, pjrho.v));
-  visc.v = vec_div(vec_mul(const_viscosity_alpha_fac.v,
-                           vec_mul(v_sig.v, vec_mul(mu_ij.v, balsara.v))),
-                   rho_ij.v);
+  visc.v = vec_div(
+      vec_mul(vec_set1(-0.25f), vec_mul(v_sig.v, vec_mul(mu_ij.v, balsara.v))),
+      rho_ij.v);
 
   /* Now, convolve with the kernel */
   visc_term.v =
@@ -937,11 +933,11 @@ runner_iact_nonsym_2_vec_force(
   rho_ij.v = vec_mul(vec_set1(0.5f), vec_add(pirho.v, pjrho.v));
   rho_ij_2.v = vec_mul(vec_set1(0.5f), vec_add(pirho.v, pjrho_2.v));
 
-  visc.v = vec_div(vec_mul(const_viscosity_alpha_fac.v,
-                           vec_mul(v_sig.v, vec_mul(mu_ij.v, balsara.v))),
-                   rho_ij.v);
+  visc.v = vec_div(
+      vec_mul(vec_set1(-0.25f), vec_mul(v_sig.v, vec_mul(mu_ij.v, balsara.v))),
+      rho_ij.v);
   visc_2.v =
-      vec_div(vec_mul(const_viscosity_alpha_fac.v,
+      vec_div(vec_mul(vec_set1(-0.25f),
                       vec_mul(v_sig_2.v, vec_mul(mu_ij_2.v, balsara_2.v))),
               rho_ij_2.v);
 
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 27aab4b897395615dcc7d0b3701e90604bd9538f..ec7d34f7ad8697f1d639ea4951011ddb06ec8833 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -200,8 +200,6 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
   io_write_attribute_s(
       h_grpsph, "Viscosity Model",
       "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
-  io_write_attribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
-  io_write_attribute_f(h_grpsph, "Viscosity beta", 3.f);
 }
 
 /**
diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index 4268a8b1a208bb1d08dfe0a838d3b59eb8cd120c..8e466daabb59482a1c2ebbaf80af30c64c4abdfe 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -451,12 +451,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
  * @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 float dt_alpha) {
+    const struct cosmology* cosmo, const struct hydro_props* hydro_props,
+    const float dt_alpha) {
 
   /* Initialise values that are used in the force loop */
   p->flux.momentum[0] = 0.0f;
diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h
index 392fb10a614f1d45227f64a459bd8d574cf6ead2..98a70aefed098243bbf2dfe08e752ee48a838d3e 100644
--- a/src/hydro/GizmoMFV/hydro.h
+++ b/src/hydro/GizmoMFV/hydro.h
@@ -476,12 +476,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
  * @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 float dt_alpha) {
+    const struct cosmology* cosmo, const struct hydro_props* hydro_props,
+    const float dt_alpha) {
 
   /* Initialise values that are used in the force loop */
   p->gravity.mflux[0] = 0.0f;
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index cfaf482278fb3624649832b2f315a039f947c7e2..1b14f05061c77681f511ff4c6ac08cc3f254b4fc 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -389,6 +389,10 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
   p->density.rho_dh = 0.f;
+  p->density.div_v = 0.f;
+  p->density.rot_v[0] = 0.f;
+  p->density.rot_v[1] = 0.f;
+  p->density.rot_v[2] = 0.f;
 }
 
 /**
@@ -424,6 +428,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   p->density.rho_dh *= h_inv_dim_plus_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 (physical) 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 (physical) velocity divergence */
+  p->density.div_v *= h_inv_dim_plus_one * a_inv2 * rho_inv;
 }
 
 /**
@@ -451,6 +466,10 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   p->density.wcount = kernel_root * h_inv_dim;
   p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
+  p->density.div_v = 0.f;
+  p->density.rot_v[0] = 0.f;
+  p->density.rot_v[1] = 0.f;
+  p->density.rot_v[2] = 0.f;
 }
 
 /**
@@ -466,12 +485,24 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @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 float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
+
+  const float fac_mu = cosmo->a_factor_mu;
+
+  /* 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->density.div_v);
 
   /* Compute the pressure */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
@@ -484,10 +515,18 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float grad_h_term =
       1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
 
+  /* Compute the Balsara switch */
+  /* Pre-multiply in the AV factor; hydro_props are not passed to the iact
+   * functions */
+  const float balsara =
+      hydro_props->viscosity.alpha * abs_div_v /
+      (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
+
   /* Update variables. */
   p->force.f = grad_h_term;
   p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
+  p->force.balsara = balsara;
 }
 
 /**
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index f716443c78ba6e7d79e09f2db345f59c6aafac65..26a854533064540630435b7a3eddbf093ae2d98a 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -87,6 +87,33 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
   pj->density.wcount += wj;
   pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+
+  /* Compute dv dot r */
+  float dv[3], curlvr[3];
+
+  const float faci = mj * wi_dx * r_inv;
+  const float facj = mi * wj_dx * r_inv;
+
+  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->density.div_v -= faci * dvdr;
+  pj->density.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];
+
+  pj->density.rot_v[0] += facj * curlvr[0];
+  pj->density.rot_v[1] += facj * curlvr[1];
+  pj->density.rot_v[2] += facj * curlvr[2];
 }
 
 /**
@@ -129,6 +156,27 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
   pi->density.wcount += wi;
   pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+  /* Compute dv dot r */
+  float dv[3], curlvr[3];
+
+  const float faci = mj * wi_dx * r_inv;
+
+  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->density.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];
 }
 
 /**
@@ -207,9 +255,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float cj = pj->force.soundspeed;
   const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
 
+  /* Grab balsara switches */
+  const float balsara_i = pi->force.balsara;
+  const float balsara_j = pj->force.balsara;
+
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+  const float visc = -0.25f * v_sig * (balsara_i + balsara_j) * mu_ij / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -330,9 +382,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float cj = pj->force.soundspeed;
   const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
 
+  /* Grab balsara switches */
+  const float balsara_i = pi->force.balsara;
+  const float balsara_j = pj->force.balsara;
+
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+  const float visc = -0.25f * v_sig * (balsara_i + balsara_j) * mu_ij / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index c33f1b9a214cf9839f1acb965b686d4a4962865c..1d14a94f2d91bf259df54c875a32bf3072ad33b6 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -124,6 +124,12 @@ struct part {
       /*! Derivative of density with respect to h */
       float rho_dh;
 
+      /*! Velocity divergence */
+      float div_v;
+
+      /*! Velocity curl */
+      float rot_v[3];
+
     } density;
 
     /**
@@ -150,6 +156,9 @@ struct part {
       /*! Time derivative of smoothing length  */
       float h_dt;
 
+      /*! Balsara switch */
+      float balsara;
+
     } force;
   };
 
diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h
index 1001dbda8b6b80a9e1c621a724350be4ec9a55c9..dee65a15758043d2cf526ea889b993c694d5dab4 100644
--- a/src/hydro/Planetary/hydro.h
+++ b/src/hydro/Planetary/hydro.h
@@ -488,14 +488,15 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @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 float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
-#ifndef PLANETARY_SPH_NO_BALSARA
   const float fac_mu = cosmo->a_factor_mu;
 
   /* Compute the norm of the curl */
@@ -505,7 +506,6 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Compute the norm of div v */
   const float abs_div_v = fabsf(p->density.div_v);
-#endif
 
   /* Compute the pressure */
   const float pressure =
@@ -527,20 +527,20 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     grad_h_term = 0.f;
   }
 
-#ifndef PLANETARY_SPH_NO_BALSARA
-  /* Compute the Balsara switch */
+    /* Compute the Balsara switch */
+#ifdef PLANETARY_SPH_NO_BALSARA
+  const float balsara = hydro_props->viscosity.alpha;
+#else
   const float balsara =
-      abs_div_v / (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
+      hydro_props->viscosity.alpha * abs_div_v /
+      (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
 #endif
 
   /* Update variables. */
   p->force.f = grad_h_term;
   p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
-
-#ifndef PLANETARY_SPH_NO_BALSARA
   p->force.balsara = balsara;
-#endif
 }
 
 /**
diff --git a/src/hydro/Planetary/hydro_iact.h b/src/hydro/Planetary/hydro_iact.h
index 149ad9dda2c2b63cc1cdca1e55280c4c7676b8b6..19ee002b85c1b0bc8ed621a029059cd02c5e670f 100644
--- a/src/hydro/Planetary/hydro_iact.h
+++ b/src/hydro/Planetary/hydro_iact.h
@@ -176,11 +176,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
                      (pi->v[1] - pj->v[1]) * dx[1] +
                      (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
 
-#ifndef PLANETARY_SPH_NO_BALSARA
   /* Balsara term */
   const float balsara_i = pi->force.balsara;
   const float balsara_j = pj->force.balsara;
-#endif
 
   /* Are the particles moving towards each other? */
   const float omega_ij = min(dvdr, 0.f);
@@ -193,12 +191,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-#ifdef PLANETARY_SPH_NO_BALSARA
-  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
-#else
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
-#endif
+  const float visc = -0.25f * 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;
@@ -300,11 +293,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
                      (pi->v[1] - pj->v[1]) * dx[1] +
                      (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
 
-#ifndef PLANETARY_SPH_NO_BALSARA
   /* Balsara term */
   const float balsara_i = pi->force.balsara;
   const float balsara_j = pj->force.balsara;
-#endif
 
   /* Are the particles moving towards each other? */
   const float omega_ij = min(dvdr, 0.f);
@@ -319,12 +310,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-#ifdef PLANETARY_SPH_NO_BALSARA
-  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
-#else
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
-#endif
+  const float visc = -0.25f * 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;
diff --git a/src/hydro/Planetary/hydro_part.h b/src/hydro/Planetary/hydro_part.h
index a40d86a53be371f3d3f50cdb9d263732cb2bdf08..4087cef62e873231a556f82869a7f6d848c8d72c 100644
--- a/src/hydro/Planetary/hydro_part.h
+++ b/src/hydro/Planetary/hydro_part.h
@@ -126,13 +126,11 @@ struct part {
       /*! Derivative of density with respect to h */
       float rho_dh;
 
-#ifndef PLANETARY_SPH_NO_BALSARA
       /*! Velocity divergence. */
       float div_v;
 
       /*! Velocity curl. */
       float rot_v[3];
-#endif
 
     } density;
 
@@ -160,10 +158,8 @@ struct part {
       /*! Time derivative of smoothing length  */
       float h_dt;
 
-#ifndef PLANETARY_SPH_NO_BALSARA
       /*! Balsara switch */
       float balsara;
-#endif
 
     } force;
   };
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index df7de2124e377ecfe41247429eba31a23575dc35..3274c1b1b58adf0fa6a4c94132fa5a87186e59ce 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -524,12 +524,14 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @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 float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
   const float fac_mu = cosmo->a_factor_mu;
 
@@ -546,7 +548,8 @@ __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 / p->h);
+      hydro_props->viscosity.alpha * abs_div_v /
+      (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu / p->h);
 
   /* Compute the "grad h" term */
   const float common_factor = p->h / (hydro_dimension * p->density.wcount);
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index 78aa4e724d5b61a125a30bd900211ea284f0db94..2ed7fe8cb8112c42de933a2ca315966e67108f0a 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -252,8 +252,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * 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;
@@ -380,8 +379,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * 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;
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 051b93e81f012106eed173679d3723dd5058827d..b16d24cfcee9407c8213b1e17465005884da6617 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -478,12 +478,14 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @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 float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
   const float fac_mu = cosmo->a_factor_mu;
 
@@ -503,7 +505,8 @@ __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 / p->h);
+      hydro_props->viscosity.alpha * abs_div_v /
+      (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu / p->h);
 
   /* Divide the pressure by the density squared to get the SPH term */
   const float rho_bar_inv = 1.f / p->rho_bar;
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
index f93e99137b8ffc47c495669047d25ab17062c593..a018b39a99be5ed691485d93bd8dfd1735378bda 100644
--- a/src/hydro/PressureEntropy/hydro_iact.h
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -263,8 +263,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
@@ -377,8 +376,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
index 3c54d0e6385480c5905aa769d51546827390951d..e9397bf6108b8bc16658157e424055274f05f23c 100644
--- a/src/hydro/PressureEntropy/hydro_io.h
+++ b/src/hydro/PressureEntropy/hydro_io.h
@@ -194,8 +194,6 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
   io_write_attribute_s(
       h_grpsph, "Viscosity Model",
       "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
-  io_write_attribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
-  io_write_attribute_f(h_grpsph, "Viscosity beta", 3.f);
 
   /* Time integration properties */
   io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 0715223850023be90294e5ccaa43c197c8eab2d4..7e38aa6b57f383564e96d9fea24730926c0ac70b 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -295,12 +295,14 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @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 float dt_alpha) {
+    const struct cosmology* cosmo, const struct hydro_props* hydro_props,
+    const float dt_alpha) {
 
   /* Initialize time step criterion variables */
   p->timestepvars.vmax = 0.0f;
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 534319b5effb6e10ce831c3504dfb0c0eb99d41f..2b1cd42055c66768e943241c75298e53e0bf75a8 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -40,6 +40,13 @@
 #define hydro_props_default_init_temp 0.f
 #define hydro_props_default_min_temp 0.f
 #define hydro_props_default_H_ionization_temperature 1e4
+#define hydro_props_default_viscosity_alpha 0.8f
+#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 */
+#define hydro_props_default_viscosity_length \
+  0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
 
 /**
  * @brief Initialize the global properties of the hydro scheme.
@@ -112,6 +119,21 @@ void hydro_props_init(struct hydro_props *p,
   p->hydrogen_mass_fraction = parser_get_opt_param_double(
       params, "SPH:H_mass_fraction", default_H_fraction);
 
+  /* Read the artificial viscosity parameters from the file, if they exist */
+  p->viscosity.alpha = parser_get_opt_param_float(
+      params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
+
+  p->viscosity.alpha_max =
+      parser_get_opt_param_float(params, "SPH:viscosity_alpha_max",
+                                 hydro_props_default_viscosity_alpha_max);
+
+  p->viscosity.alpha_min =
+      parser_get_opt_param_float(params, "SPH:viscosity_alpha_min",
+                                 hydro_props_default_viscosity_alpha_min);
+
+  p->viscosity.length = parser_get_opt_param_float(
+      params, "SPH:viscosity_length", hydro_props_default_viscosity_length);
+
   /* 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;
@@ -164,6 +186,12 @@ void hydro_props_print(const struct hydro_props *p) {
 
   message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
 
+  message(
+      "Artificial viscosity parameters set to alpha: %.3f, max: %.3f, "
+      "min: %.3f, length: %.3f.",
+      p->viscosity.alpha, p->viscosity.alpha_max, p->viscosity.alpha_min,
+      p->viscosity.length);
+
   message(
       "Hydrodynamic integration: Max change of volume: %.2f "
       "(max|dlog(h)/dt|=%f).",
@@ -221,10 +249,54 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
                        p->hydrogen_mass_fraction);
   io_write_attribute_f(h_grpsph, "Hydrogen ionization transition temperature",
                        p->hydrogen_ionization_temperature);
-  io_write_attribute_f(h_grpsph, "Alpha viscosity", const_viscosity_alpha);
+  io_write_attribute_f(h_grpsph, "Alpha viscosity", p->viscosity.alpha);
+  io_write_attribute_f(h_grpsph, "Alpha viscosity (max)",
+                       p->viscosity.alpha_max);
+  io_write_attribute_f(h_grpsph, "Alpha viscosity (min)",
+                       p->viscosity.alpha_min);
+  io_write_attribute_f(h_grpsph, "Viscosity decay length", p->viscosity.length);
+  io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
 }
 #endif
 
+/**
+ * @brief Initialises a hydro_props struct with somewhat useful values for
+ *        the automated test suite. This is not intended for production use,
+ *        but rather to fill for the purposes of mocking.
+ *
+ * @param p the struct
+ */
+void hydro_props_init_no_hydro(struct hydro_props *p) {
+  p->eta_neighbours = 1.2348;
+  p->h_tolerance = hydro_props_default_h_tolerance;
+  p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
+  const float delta_eta = p->eta_neighbours * (1.f + p->h_tolerance);
+  p->delta_neighbours =
+      (pow_dimension(delta_eta) - pow_dimension(p->eta_neighbours)) *
+      kernel_norm;
+  p->h_max = hydro_props_default_h_max;
+  p->max_smoothing_iterations = hydro_props_default_max_iterations;
+  p->CFL_condition = 0.1;
+  p->log_max_h_change = logf(powf(1.4, hydro_dimension_inv));
+
+  /* These values are inconsistent and in a production run would probably lead
+     to a crash. Again, this function is intended for mocking use in unit tests
+     and is _not_ to be used otherwise! */
+  p->minimal_temperature = hydro_props_default_min_temp;
+  p->minimal_internal_energy = 0.f;
+  p->initial_temperature = hydro_props_default_init_temp;
+  p->initial_internal_energy = 0.f;
+
+  p->hydrogen_mass_fraction = 0.755;
+  p->hydrogen_ionization_temperature =
+      hydro_props_default_H_ionization_temperature;
+
+  p->viscosity.alpha = hydro_props_default_viscosity_alpha;
+  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;
+}
+
 /**
  * @brief Write a hydro_props struct to the given FILE as a stream of bytes.
  *
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index cb3fba9724b6c2e01a616505345a99c2259a67a2..41686a67832e85c823d7aa30a30ffa6338108ee5 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -78,11 +78,27 @@ struct hydro_props {
   /*! Initial internal energy per unit mass */
   float initial_internal_energy;
 
-  /*! Primoridal hydrogen mass fraction for initial energy conversion */
+  /*! Primordial hydrogen mass fraction for initial energy conversion */
   float hydrogen_mass_fraction;
 
   /*! Temperature of the neutral to ionized transition of Hydrogen */
   float hydrogen_ionization_temperature;
+
+  /*! Artificial viscosity parameters */
+  struct {
+    /*! For the fixed, simple case. Also used to set the initial AV
+        coefficient for variable schemes. */
+    float alpha;
+
+    /*! Artificial viscosity (max) for the variable case (e.g. M&M) */
+    float alpha_max;
+
+    /*! Artificial viscosity (min) for the variable case (e.g. M&M) */
+    float alpha_min;
+
+    /*! The decay length of the artificial viscosity (used in M&M, etc.) */
+    float length;
+  } viscosity;
 };
 
 void hydro_props_print(const struct hydro_props *p);
@@ -99,4 +115,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p);
 void hydro_props_struct_dump(const struct hydro_props *p, FILE *stream);
 void hydro_props_struct_restore(const struct hydro_props *p, FILE *stream);
 
+/* Setup for tests */
+void hydro_props_init_no_hydro(struct hydro_props *p);
+
 #endif /* SWIFT_HYDRO_PROPERTIES */
diff --git a/src/map.c b/src/map.c
index 6b693885b08e1c103faa0b1b33d53cdc996c8877..68c3618fcdb10a618a97e5d1a2565d58db677cdb 100644
--- a/src/map.c
+++ b/src/map.c
@@ -191,6 +191,13 @@ void map_h_max(struct part *p, struct cell *c, void *data) {
   if (p->h > (*p2)->h) *p2 = p;
 }
 
+void map_stars_h_max(struct spart *p, struct cell *c, void *data) {
+
+  struct spart **p2 = (struct spart **)data;
+
+  if (p->h > (*p2)->h) *p2 = p;
+}
+
 /**
  * @brief Mapping function for neighbour count.
  */
diff --git a/src/map.h b/src/map.h
index 950a5fd96ebdc7177b41912b1565163f33de8701..6ad05e30df0644e1ee37b1b912bc11681ccf837c 100644
--- a/src/map.h
+++ b/src/map.h
@@ -34,6 +34,7 @@ void map_wcount_min(struct part *p, struct cell *c, void *data);
 void map_wcount_max(struct part *p, struct cell *c, void *data);
 void map_h_min(struct part *p, struct cell *c, void *data);
 void map_h_max(struct part *p, struct cell *c, void *data);
+void map_stars_h_max(struct spart *p, struct cell *c, void *data);
 void map_icount(struct part *p, struct cell *c, void *data);
 void map_dump(struct part *p, struct cell *c, void *data);
 
diff --git a/src/outputlist.c b/src/outputlist.c
index fd33370ca45f25c17ecd2cc8df622138842507f3..2ab904d4fd0b7008b324f3c37a5cab6c6b337520 100644
--- a/src/outputlist.c
+++ b/src/outputlist.c
@@ -112,6 +112,9 @@ void output_list_read_file(struct output_list *outputlist, const char *filename,
     ind += 1;
   }
 
+  /* Cleanup */
+  free(line);
+
   if (ind != outputlist->size)
     error("Did not read the correct number of output times.");
 
@@ -266,8 +269,12 @@ void output_list_print(const struct output_list *outputlist) {
 /**
  * @brief Clean an #output_list
  */
-void output_list_clean(struct output_list *outputlist) {
-  free(outputlist->times);
+void output_list_clean(struct output_list **outputlist) {
+  if (*outputlist) {
+    free((*outputlist)->times);
+    free(*outputlist);
+    *outputlist = NULL;
+  }
 }
 
 /**
diff --git a/src/outputlist.h b/src/outputlist.h
index 6045d75ea29f0aab44252835147502f3df0de20c..b7b12ca32f469c70f716553b30a15f48198f8e5e 100644
--- a/src/outputlist.h
+++ b/src/outputlist.h
@@ -58,7 +58,7 @@ void output_list_read_next_time(struct output_list *t, const struct engine *e,
 void output_list_init(struct output_list **list, const struct engine *e,
                       const char *name, double *delta_time, double *time_first);
 void output_list_print(const struct output_list *outputlist);
-void output_list_clean(struct output_list *outputlist);
+void output_list_clean(struct output_list **outputlist);
 void output_list_struct_dump(struct output_list *list, FILE *stream);
 void output_list_struct_restore(struct output_list *list, FILE *stream);
 
diff --git a/src/runner.c b/src/runner.c
index 4b19772b4b500bab49e53eda82e7cb93729f9d63..73c1e13a6cd8554ac514e37edf86f108e242e99e 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -283,7 +283,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
         for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
 
           /* Run through this cell's density interactions. */
-          for (struct link *l = finger->hydro.density; l != NULL; l = l->next) {
+          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
 
 #ifdef SWIFT_DEBUG_CHECKS
             if (l->t->ti_run < r->e->ti_current)
@@ -905,6 +905,7 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
   const int with_cosmology = (e->policy & engine_policy_cosmology);
   const double time_base = e->time_base;
   const struct cosmology *cosmo = e->cosmology;
+  const struct hydro_props *hydro_props = e->hydro_properties;
 
   TIMER_TIC;
 
@@ -943,7 +944,7 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
         }
 
         /* Compute variables required for the force loop */
-        hydro_prepare_force(p, xp, cosmo, dt_alpha);
+        hydro_prepare_force(p, xp, cosmo, hydro_props, dt_alpha);
 
         /* The particle force values are now set.  Do _NOT_
            try to read any particle density variables! */
@@ -1076,6 +1077,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
             hydro_reset_gradient(p);
 
 #else
+            const struct hydro_props *hydro_props = e->hydro_properties;
+
             /* Calculate the time-step for passing to hydro_prepare_force, used
              * for the evolution of alpha factors (i.e. those involved in the
              * artificial viscosity and thermal conduction terms) */
@@ -1095,7 +1098,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
             /* As of here, particle force variables will be set. */
 
             /* Compute variables required for the force loop */
-            hydro_prepare_force(p, xp, cosmo, dt_alpha);
+            hydro_prepare_force(p, xp, cosmo, hydro_props, dt_alpha);
 
             /* The particle force values are now set.  Do _NOT_
                try to read any particle density variables! */
@@ -1105,6 +1108,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
 #endif /* EXTRA_HYDRO_LOOP */
 
+            /* Ok, we are done with this particle */
             continue;
           }
 
@@ -1174,6 +1178,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
         hydro_reset_gradient(p);
 
 #else
+        const struct hydro_props *hydro_props = e->hydro_properties;
+
         /* Calculate the time-step for passing to hydro_prepare_force, used for
          * the evolution of alpha factors (i.e. those involved in the artificial
          * viscosity and thermal conduction terms) */
@@ -1192,7 +1198,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
         /* As of here, particle force variables will be set. */
 
         /* Compute variables required for the force loop */
-        hydro_prepare_force(p, xp, cosmo, dt_alpha);
+        hydro_prepare_force(p, xp, cosmo, hydro_props, dt_alpha);
 
         /* The particle force values are now set.  Do _NOT_
            try to read any particle density variables! */
@@ -2373,6 +2379,8 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
   const size_t nr_sparts = c->stars.count;
   const integertime_t ti_current = r->e->ti_current;
 
+  error("Need to add h_max computation");
+
   TIMER_TIC;
 
   integertime_t ti_gravity_end_min = max_nr_timesteps;
diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
index 5e656e1420134026ea6bc1d5c9c3df27887a4797..e696e4fd10853536008d9d9fafc90e6475fd291a 100644
--- a/src/runner_doiact_stars.h
+++ b/src/runner_doiact_stars.h
@@ -963,7 +963,7 @@ void runner_doself_branch_stars_density(struct runner *r, struct cell *c) {
   if (!cell_is_active_stars(c, e)) return;
 
   /* Did we mess up the recursion? */
-  if (c->hydro.h_max_old * kernel_gamma > c->dmin)
+  if (c->stars.h_max_old * kernel_gamma > c->dmin)
     error("Cell smaller than smoothing length");
 
   runner_doself_stars_density(r, c, 1);
diff --git a/src/scheduler.c b/src/scheduler.c
index 0bddcbe2c5a7d7c091f523d829b2f041725e2330..4c1c46a2c8cc6fa0697ac557e7560a1bcb128d2c 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -380,10 +380,13 @@ void scheduler_write_dependencies(struct scheduler *s, int verbose) {
     }
   }
 
-  /* Be clean */
+  /* Close the file */
   fprintf(f, "}");
   fclose(f);
+
+  /* Be clean */
   free(table);
+  free(count_rel);
 
   if (verbose)
     message("Printing task graph took %.3f %s.",
@@ -873,8 +876,6 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
  */
 static void scheduler_splittask_stars(struct task *t, struct scheduler *s) {
 
-  // LOIC: This is un-tested. Need to verify that it works.
-
   /* Iterate on this task until we're done with it. */
   int redo = 1;
   while (redo) {
@@ -1826,6 +1827,7 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
     const float gcount_i = (t->ci != NULL) ? t->ci->grav.count : 0.f;
     const float gcount_j = (t->cj != NULL) ? t->cj->grav.count : 0.f;
     const float scount_i = (t->ci != NULL) ? t->ci->stars.count : 0.f;
+    const float scount_j = (t->cj != NULL) ? t->cj->stars.count : 0.f;
 
     switch (t->type) {
       case task_type_sort:
@@ -1834,22 +1836,30 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
         break;
 
       case task_type_self:
-        // LOIC: Need to do something for stars here
         if (t->subtype == task_subtype_grav)
           cost = 1.f * (wscale * gcount_i) * gcount_i;
         else if (t->subtype == task_subtype_external_grav)
           cost = 1.f * wscale * gcount_i;
+        else if (t->subtype == task_subtype_stars_density)
+          cost = 1.f * wscale * scount_i * count_i;
         else
           cost = 1.f * (wscale * count_i) * count_i;
         break;
 
       case task_type_pair:
-        // LOIC: Need to do something for stars here
         if (t->subtype == task_subtype_grav) {
           if (t->ci->nodeID != nodeID || t->cj->nodeID != nodeID)
             cost = 3.f * (wscale * gcount_i) * gcount_j;
           else
             cost = 2.f * (wscale * gcount_i) * gcount_j;
+        } else if (t->subtype == task_subtype_stars_density) {
+          if (t->ci->nodeID != nodeID)
+            cost = 3.f * wscale * count_i * scount_j * sid_scale[t->flags];
+          else if (t->cj->nodeID != nodeID)
+            cost = 3.f * wscale * scount_i * count_j * sid_scale[t->flags];
+          else
+            cost = 2.f * wscale * (scount_i * count_j + scount_j * count_i) *
+                   sid_scale[t->flags];
         } else {
           if (t->ci->nodeID != nodeID || t->cj->nodeID != nodeID)
             cost = 3.f * (wscale * count_i) * count_j * sid_scale[t->flags];
@@ -1859,23 +1869,34 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
         break;
 
       case task_type_sub_pair:
-        // LOIC: Need to do something for stars here
-        if (t->ci->nodeID != nodeID || t->cj->nodeID != nodeID) {
 #ifdef SWIFT_DEBUG_CHECKS
-          if (t->flags < 0) error("Negative flag value!");
+        if (t->flags < 0) error("Negative flag value!");
 #endif
-          cost = 3.f * (wscale * count_i) * count_j * sid_scale[t->flags];
+        if (t->subtype == task_subtype_stars_density) {
+          if (t->ci->nodeID != nodeID) {
+            cost = 3.f * (wscale * count_i) * scount_j * sid_scale[t->flags];
+          } else if (t->cj->nodeID != nodeID) {
+            cost = 3.f * (wscale * scount_i) * count_j * sid_scale[t->flags];
+          } else {
+            cost = 2.f * wscale * (scount_i * count_j + scount_j * count_i) *
+                   sid_scale[t->flags];
+          }
+
         } else {
-#ifdef SWIFT_DEBUG_CHECKS
-          if (t->flags < 0) error("Negative flag value!");
-#endif
-          cost = 2.f * (wscale * count_i) * count_j * sid_scale[t->flags];
+          if (t->ci->nodeID != nodeID || t->cj->nodeID != nodeID) {
+            cost = 3.f * (wscale * count_i) * count_j * sid_scale[t->flags];
+          } else {
+            cost = 2.f * (wscale * count_i) * count_j * sid_scale[t->flags];
+          }
         }
         break;
 
       case task_type_sub_self:
-        // LOIC: Need to do something for stars here
-        cost = 1.f * (wscale * count_i) * count_i;
+        if (t->subtype == task_subtype_stars_density) {
+          cost = 1.f * (wscale * scount_i) * count_i;
+        } else {
+          cost = 1.f * (wscale * count_i) * count_i;
+        }
         break;
       case task_type_ghost:
         if (t->ci == t->ci->hydro.super) cost = wscale * count_i;
@@ -2523,3 +2544,61 @@ void scheduler_free_tasks(struct scheduler *s) {
   }
   s->size = 0;
 }
+
+/**
+ * @brief write down each task level
+ */
+void scheduler_write_task_level(const struct scheduler *s) {
+  /* init */
+  const int max_depth = 30;
+  const struct task *tasks = s->tasks;
+  int nr_tasks = s->nr_tasks;
+
+  /* Init counter */
+  int size = task_type_count * task_subtype_count * max_depth;
+  int *count = (int *)malloc(size * sizeof(int));
+  if (count == NULL) error("Failed to allocate memory");
+
+  for (int i = 0; i < size; i++) count[i] = 0;
+
+  /* Count tasks */
+  for (int i = 0; i < nr_tasks; i++) {
+    const struct task *t = &tasks[i];
+    if (t->ci) {
+
+      if ((int)t->ci->depth >= max_depth)
+        error("Cell is too deep, you need to increase max_depth");
+
+      int ind = t->type * task_subtype_count * max_depth;
+      ind += t->subtype * max_depth;
+      ind += (int)t->ci->depth;
+
+      count[ind] += 1;
+    }
+  }
+
+  /* Open file */
+  char filename[200] = "task_level.txt";
+  FILE *f = fopen(filename, "w");
+  if (f == NULL) error("Error opening task level file.");
+
+  /* Print header */
+  fprintf(f, "# task_type, task_subtype, depth, count\n");
+
+  /* Print tasks level */
+  for (int i = 0; i < size; i++) {
+    if (count[i] == 0) continue;
+
+    int type = i / (task_subtype_count * max_depth);
+    int subtype = i - task_subtype_count * max_depth * type;
+    subtype /= max_depth;
+    int depth = i - task_subtype_count * max_depth * type;
+    depth -= subtype * max_depth;
+    fprintf(f, "%s %s %i %i\n", taskID_names[type], subtaskID_names[subtype],
+            depth, count[i]);
+  }
+
+  /* clean up */
+  fclose(f);
+  free(count);
+}
diff --git a/src/scheduler.h b/src/scheduler.h
index 1a75544de12b8402e553e3ae2b84e2d8a65c56e8..f1e130c6ce2a8538b0126e86ee0cbd79cf5a0e0d 100644
--- a/src/scheduler.h
+++ b/src/scheduler.h
@@ -173,5 +173,6 @@ void scheduler_print_tasks(const struct scheduler *s, const char *fileName);
 void scheduler_clean(struct scheduler *s);
 void scheduler_free_tasks(struct scheduler *s);
 void scheduler_write_dependencies(struct scheduler *s, int verbose);
+void scheduler_write_task_level(const struct scheduler *s);
 
 #endif /* SWIFT_SCHEDULER_H */
diff --git a/src/space.c b/src/space.c
index da81a667fa64782106aea03526410bde1cb03f6e..5ee6e4c3c1934e6c3c6c4995e59ad79140c4139d 100644
--- a/src/space.c
+++ b/src/space.c
@@ -176,6 +176,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->grav.mm = NULL;
     c->hydro.dx_max_part = 0.0f;
     c->hydro.dx_max_sort = 0.0f;
+    c->stars.dx_max_part = 0.f;
     c->hydro.sorted = 0;
     c->hydro.count = 0;
     c->hydro.updated = 0;
@@ -274,6 +275,7 @@ void space_free_cells(struct space *s) {
 void space_regrid(struct space *s, int verbose) {
 
   const size_t nr_parts = s->nr_parts;
+  const size_t nr_sparts = s->nr_sparts;
   const ticks tic = getticks();
   const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
 
@@ -287,6 +289,9 @@ void space_regrid(struct space *s, int verbose) {
         if (c->hydro.h_max > h_max) {
           h_max = s->cells_top[k].hydro.h_max;
         }
+        if (c->stars.h_max > h_max) {
+          h_max = s->cells_top[k].stars.h_max;
+        }
       }
     } else if (s->cells_top != NULL) {
       for (int k = 0; k < s->nr_cells; k++) {
@@ -294,11 +299,17 @@ void space_regrid(struct space *s, int verbose) {
         if (c->nodeID == engine_rank && c->hydro.h_max > h_max) {
           h_max = s->cells_top[k].hydro.h_max;
         }
+        if (c->nodeID == engine_rank && c->stars.h_max > h_max) {
+          h_max = s->cells_top[k].stars.h_max;
+        }
       }
     } else {
       for (size_t k = 0; k < nr_parts; k++) {
         if (s->parts[k].h > h_max) h_max = s->parts[k].h;
       }
+      for (size_t k = 0; k < nr_sparts; k++) {
+        if (s->sparts[k].h > h_max) h_max = s->sparts[k].h;
+      }
     }
   }
 
@@ -1911,6 +1922,7 @@ void space_split_recursive(struct space *s, struct cell *c,
   const int depth = c->depth;
   int maxdepth = 0;
   float h_max = 0.0f;
+  float stars_h_max = 0.f;
   integertime_t ti_hydro_end_min = max_nr_timesteps, ti_hydro_end_max = 0,
                 ti_hydro_beg_max = 0;
   integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_end_max = 0,
@@ -2013,6 +2025,8 @@ void space_split_recursive(struct space *s, struct cell *c,
       cp->hydro.h_max = 0.f;
       cp->hydro.dx_max_part = 0.f;
       cp->hydro.dx_max_sort = 0.f;
+      cp->stars.h_max = 0.f;
+      cp->stars.dx_max_part = 0.f;
       cp->nodeID = c->nodeID;
       cp->parent = c;
       cp->super = NULL;
@@ -2061,6 +2075,7 @@ void space_split_recursive(struct space *s, struct cell *c,
 
         /* Update the cell-wide properties */
         h_max = max(h_max, cp->hydro.h_max);
+        stars_h_max = max(h_max, cp->stars.h_max);
         ti_hydro_end_min = min(ti_hydro_end_min, cp->hydro.ti_end_min);
         ti_hydro_end_max = max(ti_hydro_end_max, cp->hydro.ti_end_max);
         ti_hydro_beg_max = max(ti_hydro_beg_max, cp->hydro.ti_beg_max);
@@ -2230,6 +2245,12 @@ void space_split_recursive(struct space *s, struct cell *c,
 #endif
       gravity_time_bin_min = min(gravity_time_bin_min, sparts[k].time_bin);
       gravity_time_bin_max = max(gravity_time_bin_max, sparts[k].time_bin);
+      stars_h_max = max(stars_h_max, sparts[k].h);
+
+      /* Reset x_diff */
+      sparts[k].x_diff[0] = 0.f;
+      sparts[k].x_diff[1] = 0.f;
+      sparts[k].x_diff[2] = 0.f;
     }
 
     /* Convert into integer times */
@@ -2278,6 +2299,7 @@ void space_split_recursive(struct space *s, struct cell *c,
   c->grav.ti_end_min = ti_gravity_end_min;
   c->grav.ti_end_max = ti_gravity_end_max;
   c->grav.ti_beg_max = ti_gravity_beg_max;
+  c->stars.h_max = stars_h_max;
   c->maxdepth = maxdepth;
 
   /* Set ownership according to the start of the parts array. */
diff --git a/src/stars/Default/stars_part.h b/src/stars/Default/stars_part.h
index fe5be3c7f6f190adea2c81ee14c0cad7a10fbd2c..1922cc0e260a3c0f2052649a87252019514e637f 100644
--- a/src/stars/Default/stars_part.h
+++ b/src/stars/Default/stars_part.h
@@ -38,6 +38,9 @@ struct spart {
   /*! Particle position. */
   double x[3];
 
+  /* Offset between current position and position at last tree rebuild. */
+  float x_diff[3];
+
   /*! Particle velocity. */
   float v[3];
 
diff --git a/src/velociraptor_dummy.c b/src/velociraptor_dummy.c
index 9b364709b8eb12ce42eead0d10dae8e87d907d49..8f14a3230d341993122f09f2bccf3d8232550fd9 100644
--- a/src/velociraptor_dummy.c
+++ b/src/velociraptor_dummy.c
@@ -24,21 +24,31 @@
 #include <stddef.h>
 
 /* Local includes. */
+#include "error.h"
 #include "swift_velociraptor_part.h"
 #include "velociraptor_interface.h"
 
 /* Dummy VELOCIraptor interface for testing compilation without linking the
  * actual VELOCIraptor library. */
 #ifdef HAVE_DUMMY_VELOCIRAPTOR
+struct cosmoinfo {};
+struct unitinfo {};
+struct cell_loc {};
+struct siminfo {};
+
 int InitVelociraptor(char *config_name, char *output_name,
                      struct cosmoinfo cosmo_info, struct unitinfo unit_info,
                      struct siminfo sim_info) {
+
+  error("This is only a dummy. Call the real one!");
   return 0;
 }
 int InvokeVelociraptor(const size_t num_gravity_parts,
                        const size_t num_hydro_parts,
                        struct swift_vel_part *swift_parts,
                        const int *cell_node_ids, char *output_name) {
+
+  error("This is only a dummy. Call the real one!");
   return 0;
 }
 #endif /* HAVE_DUMMY_VELOCIRAPTOR */
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 0b52bf68734e398db515e204fa3cec60d5c6f314..93a85bea87eb1b7c204f0cf9e6ea37ecea6d18f4 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -31,19 +31,12 @@
 #include "swift.h"
 
 #if defined(WITH_VECTORIZATION)
-#define DOSELF2 runner_doself2_force_vec
-#define DOPAIR2 runner_dopair2_branch_force
 #define DOSELF2_NAME "runner_doself2_force_vec"
 #define DOPAIR2_NAME "runner_dopair2_force_vec"
 #endif
 
-#ifndef DOSELF2
-#define DOSELF2 runner_doself2_force
+#ifndef DOSELF2_NAME
 #define DOSELF2_NAME "runner_doself2_density"
-#endif
-
-#ifndef DOPAIR2
-#define DOPAIR2 runner_dopair2_branch_force
 #define DOPAIR2_NAME "runner_dopair2_force"
 #endif
 
@@ -590,6 +583,7 @@ int main(int argc, char *argv[]) {
   prog_const.const_newton_G = 1.f;
 
   struct hydro_props hp;
+  hydro_props_init_no_hydro(&hp);
   hp.eta_neighbours = h;
   hp.h_tolerance = 1e0;
   hp.h_max = FLT_MAX;
@@ -669,14 +663,13 @@ int main(int argc, char *argv[]) {
     for (int j = 0; j < 125; ++j)
       runner_do_sort(&runner, cells[j], 0x1FFF, 0, 0);
 
-/* Do the density calculation */
-#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
+      /* Do the density calculation */
 
 /* Initialise the particle cache. */
 #ifdef WITH_VECTORIZATION
     runner.ci_cache.count = 0;
-    cache_init(&runner.ci_cache, 512);
     runner.cj_cache.count = 0;
+    cache_init(&runner.ci_cache, 512);
     cache_init(&runner.cj_cache, 512);
 #endif
 
@@ -714,18 +707,15 @@ int main(int argc, char *argv[]) {
     for (int j = 0; j < 27; ++j)
       runner_doself1_density(&runner, inner_cells[j]);
 
-#endif
-
     /* Ghost to finish everything on the central cells */
     for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j], 0);
 
-/* Do the force calculation */
-#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
+      /* Do the force calculation */
 
 #ifdef WITH_VECTORIZATION
     /* Initialise the cache. */
-    runner.ci_cache.count = 0;
-    runner.cj_cache.count = 0;
+    cache_clean(&runner.ci_cache);
+    cache_clean(&runner.cj_cache);
     cache_init(&runner.ci_cache, 512);
     cache_init(&runner.cj_cache, 512);
 #endif
@@ -742,7 +732,7 @@ int main(int argc, char *argv[]) {
 
             const ticks sub_tic = getticks();
 
-            DOPAIR2(&runner, main_cell, cj);
+            runner_dopair2_branch_force(&runner, main_cell, cj);
 
             timings[ctr++] += getticks() - sub_tic;
           }
@@ -753,10 +743,9 @@ int main(int argc, char *argv[]) {
     ticks self_tic = getticks();
 
     /* And now the self-interaction for the main cell */
-    DOSELF2(&runner, main_cell);
+    runner_doself2_force(&runner, main_cell);
 
     timings[26] += getticks() - self_tic;
-#endif
 
     /* Finally, give a gentle kick */
     runner_do_end_force(&runner, main_cell, 0);
@@ -802,18 +791,17 @@ int main(int argc, char *argv[]) {
 
   const ticks tic = getticks();
 
-/* Kick the central cell */
-// runner_do_kick1(&runner, main_cell, 0);
+  /* Kick the central cell */
+  // runner_do_kick1(&runner, main_cell, 0);
 
-/* And drift it */
-// runner_do_drift_particles(&runner, main_cell, 0);
+  /* And drift it */
+  // runner_do_drift_particles(&runner, main_cell, 0);
 
-/* Initialise the particles */
-// for (int j = 0; j < 125; ++j) runner_do_drift_particles(&runner, cells[j],
-// 0);
+  /* Initialise the particles */
+  // for (int j = 0; j < 125; ++j) runner_do_drift_particles(&runner, cells[j],
+  // 0);
 
-/* Do the density calculation */
-#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
+  /* Do the density calculation */
 
   /* Run all the pairs (only once !)*/
   for (int i = 0; i < 5; i++) {
@@ -848,13 +836,10 @@ int main(int argc, char *argv[]) {
   /* And now the self-interaction for the central cells*/
   for (int j = 0; j < 27; ++j) self_all_density(&runner, inner_cells[j]);
 
-#endif
-
   /* Ghost to finish everything on the central cells */
   for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j], 0);
 
-/* Do the force calculation */
-#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
+  /* Do the force calculation */
 
   /* Do the pairs (for the central 27 cells) */
   for (int i = 1; i < 4; i++) {
@@ -871,8 +856,6 @@ int main(int argc, char *argv[]) {
   /* And now the self-interaction for the main cell */
   self_all_force(&runner, main_cell);
 
-#endif
-
   /* Finally, give a gentle kick */
   runner_do_end_force(&runner, main_cell, 0);
   // runner_do_kick2(&runner, main_cell, 0);
@@ -890,5 +873,10 @@ int main(int argc, char *argv[]) {
   for (int i = 0; i < 125; ++i) clean_up(cells[i]);
   free(solution);
 
+#ifdef WITH_VECTORIZATION
+  cache_clean(&runner.ci_cache);
+  cache_clean(&runner.cj_cache);
+#endif
+
   return 0;
 }
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 0ba724c7da8051c10bb8205f2d9ad2ecbf3753c3..ba98d91a4250865f301bd96594b1364d95034bbb 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -456,6 +456,7 @@ int main(int argc, char *argv[]) {
   space.dim[2] = 3.;
 
   struct hydro_props hp;
+  hydro_props_init_no_hydro(&hp);
   hp.eta_neighbours = h;
   hp.h_tolerance = 1e0;
   hp.h_max = FLT_MAX;
@@ -615,5 +616,10 @@ int main(int argc, char *argv[]) {
   /* Clean things to make the sanitizer happy ... */
   for (int i = 0; i < 27; ++i) clean_up(cells[i]);
 
+#ifdef WITH_VECTORIZATION
+  cache_clean(&runner.ci_cache);
+  cache_clean(&runner.cj_cache);
+#endif
+
   return 0;
 }
diff --git a/tests/testActivePair.c b/tests/testActivePair.c
index a96e4f661c06b75c14fd3b96fc3d28376b0b959e..dd7b4c6277142f4e074fc154aac455c355200816 100644
--- a/tests/testActivePair.c
+++ b/tests/testActivePair.c
@@ -33,7 +33,8 @@
 
 /* Typdef function pointer for interaction function. */
 typedef void (*interaction_func)(struct runner *, struct cell *, struct cell *);
-typedef void (*init_func)(struct cell *, const struct cosmology *);
+typedef void (*init_func)(struct cell *, const struct cosmology *,
+                          const struct hydro_props *);
 typedef void (*finalise_func)(struct cell *, const struct cosmology *);
 
 /**
@@ -169,8 +170,8 @@ void clean_up(struct cell *ci) {
 /**
  * @brief Initializes all particles field to be ready for a density calculation
  */
-void zero_particle_fields_density(struct cell *c,
-                                  const struct cosmology *cosmo) {
+void zero_particle_fields_density(struct cell *c, const struct cosmology *cosmo,
+                                  const struct hydro_props *hydro_props) {
   for (int pid = 0; pid < c->hydro.count; pid++) {
     hydro_init_part(&c->hydro.parts[pid], NULL);
   }
@@ -179,7 +180,8 @@ void zero_particle_fields_density(struct cell *c,
 /**
  * @brief Initializes all particles field to be ready for a force calculation
  */
-void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo) {
+void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo,
+                                const struct hydro_props *hydro_props) {
   for (int pid = 0; pid < c->hydro.count; pid++) {
     struct part *p = &c->hydro.parts[pid];
     struct xpart *xp = &c->hydro.xparts[pid];
@@ -219,7 +221,7 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo) {
 #endif /* PRESSURE-ENERGY */
 
     /* And prepare for a round of force tasks. */
-    hydro_prepare_force(p, xp, cosmo, 0.);
+    hydro_prepare_force(p, xp, cosmo, hydro_props, 0.);
     hydro_reset_acceleration(p);
   }
 }
@@ -299,8 +301,8 @@ void test_pair_interactions(struct runner *runner, struct cell **ci,
   runner_do_sort(runner, *cj, 0x1FFF, 0, 0);
 
   /* Zero the fields */
-  init(*ci, runner->e->cosmology);
-  init(*cj, runner->e->cosmology);
+  init(*ci, runner->e->cosmology, runner->e->hydro_properties);
+  init(*cj, runner->e->cosmology, runner->e->hydro_properties);
 
   /* Run the test */
   vec_interaction(runner, *ci, *cj);
@@ -315,8 +317,8 @@ void test_pair_interactions(struct runner *runner, struct cell **ci,
   /* Now perform a brute-force version for accuracy tests */
 
   /* Zero the fields */
-  init(*ci, runner->e->cosmology);
-  init(*cj, runner->e->cosmology);
+  init(*ci, runner->e->cosmology, runner->e->hydro_properties);
+  init(*cj, runner->e->cosmology, runner->e->hydro_properties);
 
   /* Run the brute-force test */
   serial_interaction(runner, *ci, *cj);
@@ -487,6 +489,7 @@ int main(int argc, char *argv[]) {
   struct space space;
   struct engine engine;
   struct cosmology cosmo;
+  struct hydro_props hydro_props;
   struct runner *runner;
   char c;
   static long long partId = 0;
@@ -571,6 +574,8 @@ int main(int argc, char *argv[]) {
 
   cosmology_init_no_cosmo(&cosmo);
   engine.cosmology = &cosmo;
+  hydro_props_init_no_hydro(&hydro_props);
+  engine.hydro_properties = &hydro_props;
 
   if (posix_memalign((void **)&runner, SWIFT_STRUCT_ALIGNMENT,
                      sizeof(struct runner)) != 0) {
diff --git a/tests/testAdiabaticIndex.c b/tests/testAdiabaticIndex.c
index 60ecefa264f48bed2d4df205766dc392a1a03d0f..6aa794207f0e23e6a26060f3ef98b7ee841d7a32 100644
--- a/tests/testAdiabaticIndex.c
+++ b/tests/testAdiabaticIndex.c
@@ -34,7 +34,8 @@
  */
 void check_value(float a, float b, const char* s) {
   if (fabsf(a - b) / fabsf(a + b) > 1.e-6f)
-    error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s);
+    error("Values are inconsistent: %12.15e %12.15e rel=%e (%s)!", a, b,
+          fabsf(a - b) / fabsf(a + b), s);
 }
 
 /**
@@ -77,36 +78,61 @@ void check_constants(void) {
 void check_functions(float x) {
 
   float val_a, val_b;
+  const double xx = x;
+
+#if defined(HYDRO_GAMMA_5_3)
+#define hydro_gamma_d (5. / 3.)
+#elif defined(HYDRO_GAMMA_7_5)
+#define hydro_gamma_d (7. / 5.)
+#elif defined(HYDRO_GAMMA_4_3)
+#define hydro_gamma_d (4. / 3.)
+#elif defined(HYDRO_GAMMA_2_1)
+#define hydro_gamma_d (2. / 1.)
+#else
+#error "Need to choose an adiabatic index!"
+#endif
+
+  val_a = pow(xx, hydro_gamma_d);
+  val_b = pow_gamma(x);
+  check_value(val_a, val_b, "x^gamma");
+
+  val_a = pow(xx, hydro_gamma_d - 1.0);
+  val_b = pow_gamma_minus_one(x);
+  check_value(val_a, val_b, "x^(gamma - 1)");
+
+  val_a = pow(xx, -(hydro_gamma_d - 1.0));
+  val_b = pow_minus_gamma_minus_one(x);
+  check_value(val_a, val_b, "x^(-(gamma - 1))");
 
-  val_a = powf(x, -hydro_gamma);
+  val_a = pow(xx, -hydro_gamma_d);
   val_b = pow_minus_gamma(x);
   check_value(val_a, val_b, "x^(-gamma)");
 
-  val_a = powf(x, 2.0f / (hydro_gamma - 1.0f));
+  val_a = pow(xx, 2.0 / (hydro_gamma_d - 1.0));
   val_b = pow_two_over_gamma_minus_one(x);
   check_value(val_a, val_b, "x^(2/(gamma-1))");
 
-  val_a = powf(x, 2.0f * hydro_gamma / (hydro_gamma - 1.0f));
+  val_a = pow(xx, 2.0 * hydro_gamma_d / (hydro_gamma_d - 1.0));
   val_b = pow_two_gamma_over_gamma_minus_one(x);
   check_value(val_a, val_b, "x^((2 gamma)/(gamma-1))");
 
-  val_a = powf(x, 0.5f * (hydro_gamma - 1.0f) / hydro_gamma);
+  val_a = pow(xx, (hydro_gamma_d - 1.0) / (2.0 * hydro_gamma_d));
   val_b = pow_gamma_minus_one_over_two_gamma(x);
   check_value(val_a, val_b, "x^((gamma-1)/(2 gamma))");
 
-  val_a = powf(x, -0.5f * (hydro_gamma + 1.0f) / hydro_gamma);
+  val_a = pow(xx, -(hydro_gamma_d + 1.0) / (2.0 * hydro_gamma_d));
   val_b = pow_minus_gamma_plus_one_over_two_gamma(x);
   check_value(val_a, val_b, "x^(-(gamma+1)/(2 gamma))");
 
-  val_a = powf(x, 1.0f / hydro_gamma);
+  val_a = pow(xx, 1.0 / hydro_gamma_d);
   val_b = pow_one_over_gamma(x);
   check_value(val_a, val_b, "x^(1/gamma)");
 
-  val_a = powf(x, 3.f * hydro_gamma - 2.f);
+  val_a = pow(xx, 3. * hydro_gamma_d - 2.);
   val_b = pow_three_gamma_minus_two(x);
   check_value(val_a, val_b, "x^(3gamma - 2)");
 
-  val_a = powf(x, (3.f * hydro_gamma - 5.f) / 2.f);
+  val_a = pow(xx, (3. * hydro_gamma_d - 5.) / 2.);
   val_b = pow_three_gamma_minus_five_over_two(x);
   check_value(val_a, val_b, "x^((3gamma - 5)/2)");
 }
diff --git a/tests/testCbrt.c b/tests/testCbrt.c
index b608f9a00d619570c298f4123038f930584a245c..3663e0e19ad2a5ad35d67703e00f5c0309a3eb00 100644
--- a/tests/testCbrt.c
+++ b/tests/testCbrt.c
@@ -125,5 +125,6 @@ int main(int argc, char *argv[]) {
   message("x * icbrtf   took %9.3f %s (acc = %18.11e).",
           clocks_from_ticks(getticks() - tic_ours), clocks_getunit(), acc_ours);
 
+  free(data);
   return 0;
 }
diff --git a/tests/testOutputList.c b/tests/testOutputList.c
index b7df197405ee095cf9bf0a63e8cf7f00585f269f..cd9c62c3abc7462406950e24fa330788b76ed249 100644
--- a/tests/testOutputList.c
+++ b/tests/testOutputList.c
@@ -76,7 +76,7 @@ void test_no_cosmo(struct engine *e, char *name, int with_assert) {
     output_time = (double)(ti_next * e->time_base) + e->time_begin;
   }
 
-  output_list_clean(list);
+  output_list_clean(&list);
 };
 
 void test_cosmo(struct engine *e, char *name, int with_assert) {
@@ -112,7 +112,7 @@ void test_cosmo(struct engine *e, char *name, int with_assert) {
     output_time = (double)exp(ti_next * e->time_base) * e->cosmology->a_begin;
   }
 
-  output_list_clean(list);
+  output_list_clean(&list);
 };
 
 int main(int argc, char *argv[]) {
@@ -151,6 +151,8 @@ int main(int argc, char *argv[]) {
   test_cosmo(&e, "ScaleFactor", with_assert);
   test_cosmo(&e, "Time", without_assert);
 
+  cosmology_clean(&cosmo);
+
   /* Write message and leave */
   message("Test done");
   return 0;
diff --git a/tests/testParser.c b/tests/testParser.c
index 3944e86fa19a1f623623383eabefe1094bf5addf..84ce70ff44fad0482573c740d5a174285655c08d 100644
--- a/tests/testParser.c
+++ b/tests/testParser.c
@@ -114,6 +114,8 @@ int main(int argc, char *argv[]) {
   int haveoptwords1 = parser_get_opt_param_string_array(
       &param_file, "Simulation:optwords", &nvar_result, &var_result, noptwords,
       optwords);
+  parser_free_param_string_array(nvar_result, var_result);
+
   /* Check if we can read it again */
   int haveoptwords2 = parser_get_opt_param_string_array(
       &param_file, "Simulation:optwords", &nvar_result, &var_result, noptwords,
diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c
index 5092df2d2a2df2d8fc6e1f3c8bbdc4432d604df9..064c86d42f8df907d1ffaaab164b6a2f8b534b19 100644
--- a/tests/testPotentialPair.c
+++ b/tests/testPotentialPair.c
@@ -439,5 +439,10 @@ int main(int argc, char *argv[]) {
   free(cj.grav.multipole);
   free(ci.grav.parts);
   free(cj.grav.parts);
+
+  /* Clean up the caches */
+  gravity_cache_clean(&r.ci_gravity_cache);
+  gravity_cache_clean(&r.cj_gravity_cache);
+
   return 0;
 }
diff --git a/tests/testPotentialSelf.c b/tests/testPotentialSelf.c
index 08a0870c61d128b62b8b72f90ecc42317acd895c..10eb499570a591daaf0de2e011f2346077905e8e 100644
--- a/tests/testPotentialSelf.c
+++ b/tests/testPotentialSelf.c
@@ -200,5 +200,10 @@ int main(int argc, char *argv[]) {
   }
 
   free(c.grav.parts);
+
+  /* Clean up the caches */
+  gravity_cache_clean(&r.ci_gravity_cache);
+
+  /* All done! */
   return 0;
 }