diff --git a/examples/CosmoVolume/run.sh b/examples/CosmoVolume/run.sh
index 412456c3cb1ae9b869afa52b1046747e32c2eefe..e128aa8f88cc30e13c420bda2088c9958c5bfc31 100755
--- a/examples/CosmoVolume/run.sh
+++ b/examples/CosmoVolume/run.sh
@@ -7,4 +7,4 @@ then
     ./getIC.sh
 fi
 
-../swift -s -t 16 cosmoVolume.yml
+../swift -s -t 16 cosmoVolume.yml 2>&1 | tee output.log
diff --git a/examples/EAGLE_12/run.sh b/examples/EAGLE_12/run.sh
index 58c92adc0a1afadf0d15f81613d749f6e736f211..21b965050ab1c786fd08c0dfd28b1612ac09e9e5 100755
--- a/examples/EAGLE_12/run.sh
+++ b/examples/EAGLE_12/run.sh
@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../swift -s -t 16 eagle_12.yml
+../swift -s -t 16 eagle_12.yml 2>&1 | tee output.log
 
diff --git a/examples/EAGLE_25/run.sh b/examples/EAGLE_25/run.sh
index f232cf98f772ab9900bd22b635090e3c103749b5..21d0c7baf4482bdc6df83c552253a02d47432ba2 100755
--- a/examples/EAGLE_25/run.sh
+++ b/examples/EAGLE_25/run.sh
@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../swift -s -t 16 eagle_25.yml
+../swift -s -t 16 eagle_25.yml 2>&1 | tee output.log
 
diff --git a/examples/EAGLE_50/run.sh b/examples/EAGLE_50/run.sh
index 87d522be6bd95fa503d437af861a68ce8bf85b4c..253e7e6eb11a04e1731c605e9ea067ba2bd13221 100755
--- a/examples/EAGLE_50/run.sh
+++ b/examples/EAGLE_50/run.sh
@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../swift -s -t 16 eagle_50.yml
+../swift -s -t 16 eagle_50.yml 2>&1 | tee output.log
 
diff --git a/examples/ExternalPointMass/run.sh b/examples/ExternalPointMass/run.sh
index 9bbe03738b472f48715ae1875dfd462ef577b3d9..9f90ca395a5c8cf83e67928b3fdbd4d8529ac254 100755
--- a/examples/ExternalPointMass/run.sh
+++ b/examples/ExternalPointMass/run.sh
@@ -7,4 +7,4 @@ then
     python makeIC.py 10000
 fi
 
-../swift -g -t 2 externalPointMass.yml
+../swift -g -t 2 externalPointMass.yml 2>&1 | tee output.log
diff --git a/examples/Glass/glass_50000.hdf5 b/examples/Glass/glass_50000.hdf5
deleted file mode 100644
index ad209b851f04ee26e527fe925b71322a2740bb55..0000000000000000000000000000000000000000
Binary files a/examples/Glass/glass_50000.hdf5 and /dev/null differ
diff --git a/examples/GreshoVortex_2D/run.sh b/examples/GreshoVortex_2D/run.sh
index 81c016f4775a1e5d693e3c2f78c7d1e5dc532701..6d537bcc96c68385434f685bd551a2d423f469e0 100755
--- a/examples/GreshoVortex_2D/run.sh
+++ b/examples/GreshoVortex_2D/run.sh
@@ -13,7 +13,7 @@ then
 fi
 
 # Run SWIFT
-../swift -s -t 1 gresho.yml
+../swift -s -t 1 gresho.yml 2>&1 | tee output.log
 
 # Plot the solution
 python plotSolution.py 11
diff --git a/examples/IsothermalPotential/GravityOnly/legend.pro b/examples/IsothermalPotential/GravityOnly/legend.pro
deleted file mode 100755
index e510bdfcd7dcabed796cb090c6fa7d54536608b6..0000000000000000000000000000000000000000
--- a/examples/IsothermalPotential/GravityOnly/legend.pro
+++ /dev/null
@@ -1,459 +0,0 @@
-;+
-; NAME:
-;       LEGEND
-; PURPOSE:
-;       Create an annotation legend for a plot.
-; EXPLANATION:
-;       This procedure makes a legend for a plot.  The legend can contain
-;       a mixture of symbols, linestyles, Hershey characters (vectorfont),
-;       and filled polygons (usersym).  A test procedure, legendtest.pro,
-;       shows legend's capabilities.  Placement of the legend is controlled
-;       with keywords like /right, /top, and /center or by using a position
-;       keyword for exact placement (position=[x,y]) or via mouse (/position).
-; CALLING SEQUENCE:
-;       LEGEND [,items][,keyword options]
-; EXAMPLES:
-;       The call:
-;               legend,['Plus sign','Asterisk','Period'],psym=[1,2,3]
-;         produces:
-;               -----------------
-;               |               |
-;               |  + Plus sign  |
-;               |  * Asterisk   |
-;               |  . Period     |
-;               |               |
-;               -----------------
-;         Each symbol is drawn with a plots command, so they look OK.
-;         Other examples are given in optional output keywords.
-;
-;       lines = indgen(6)                       ; for line styles
-;       items = 'linestyle '+strtrim(lines,2)   ; annotations
-;       legend,items,linestyle=lines            ; vertical legend---upper left
-;       items = ['Plus sign','Asterisk','Period']
-;       sym = [1,2,3]
-;       legend,items,psym=sym                   ; ditto except using symbols
-;       legend,items,psym=sym,/horizontal       ; horizontal format
-;       legend,items,psym=sym,box=0             ; sans border
-;       legend,items,psym=sym,delimiter='='     ; embed '=' betw psym & text
-;       legend,items,psym=sym,margin=2          ; 2-character margin
-;       legend,items,psym=sym,position=[x,y]    ; upper left in data coords
-;       legend,items,psym=sym,pos=[x,y],/norm   ; upper left in normal coords
-;       legend,items,psym=sym,pos=[x,y],/device ; upper left in device coords
-;       legend,items,psym=sym,/position         ; interactive position
-;       legend,items,psym=sym,/right            ; at upper right
-;       legend,items,psym=sym,/bottom           ; at lower left
-;       legend,items,psym=sym,/center           ; approximately near center
-;       legend,items,psym=sym,number=2          ; plot two symbols, not one
-;       legend,items,/fill,psym=[8,8,8],colors=[10,20,30]; 3 filled squares
-; INPUTS:
-;       items = text for the items in the legend, a string array.
-;               For example, items = ['diamond','asterisk','square'].
-;               You can omit items if you don't want any text labels.
-; OPTIONAL INPUT KEYWORDS:
-;
-;       linestyle = array of linestyle numbers  If linestyle[i] < 0, then omit
-;               ith symbol or line to allow a multi-line entry.     If 
-;               linestyle = -99 then text will be left-justified.  
-;       psym = array of plot symbol numbers.  If psym[i] is negative, then a
-;               line connects pts for ith item.  If psym[i] = 8, then the
-;               procedure usersym is called with vertices define in the
-;               keyword usersym.   If psym[i] = 88, then use the previously
-;               defined user symbol
-;       vectorfont = vector-drawn characters for the sym/line column, e.g.,
-;               ['!9B!3','!9C!3','!9D!3'] produces an open square, a checkmark,
-;               and a partial derivative, which might have accompanying items
-;               ['BOX','CHECK','PARTIAL DERIVATIVE'].
-;               There is no check that !p.font is set properly, e.g., -1 for
-;               X and 0 for PostScript.  This can produce an error, e.g., use
-;               !20 with PostScript and !p.font=0, but allows use of Hershey
-;               *AND* PostScript fonts together.
-;       N. B.: Choose any of linestyle, psym, and/or vectorfont.  If none is
-;               present, only the text is output.  If more than one
-;               is present, all need the same number of elements, and normal
-;               plot behaviour occurs.
-;               By default, if psym is positive, you get one point so there is
-;               no connecting line.  If vectorfont[i] = '',
-;               then plots is called to make a symbol or a line, but if
-;               vectorfont[i] is a non-null string, then xyouts is called.
-;       /help = flag to print header
-;       /horizontal = flag to make the legend horizontal
-;       /vertical = flag to make the legend vertical (D=vertical)
-;       box = flag to include/omit box around the legend (D=include)
-;       clear = flag to clear the box area before drawing the legend
-;       delimiter = embedded character(s) between symbol and text (D=none)
-;       colors = array of colors for plot symbols/lines (D=!P.color)
-;       textcolors = array of colors for text (D=!P.color)
-;       margin = margin around text measured in characters and lines
-;       spacing = line spacing (D=bit more than character height)
-;       pspacing = psym spacing (D=3 characters) (when number of symbols is
-;             greater than 1)
-;       charsize = just like !p.charsize for plot labels
-;       charthick = just like !p.charthick for plot labels
-;       thick = array of line thickness numbers (D = !P.thick), if used, then 
-;               linestyle must also be specified
-;       position = data coordinates of the /top (D) /left (D) of the legend
-;       normal = use normal coordinates for position, not data
-;       device = use device coordinates for position, not data
-;       number = number of plot symbols to plot or length of line (D=1)
-;       usersym = 2-D array of vertices, cf. usersym in IDL manual. 
-;             (/USERSYM =square, default is to use existing USERSYM definition)
-;       /fill = flag to fill the usersym
-;       /left_legend = flag to place legend snug against left side of plot
-;                 window (D)
-;       /right_legend = flag to place legend snug against right side of plot
-;               window.    If /right,pos=[x,y], then x is position of RHS and
-;               text runs right-to-left.
-;       /top_legend = flag to place legend snug against top of plot window (D)
-;       /bottom = flag to place legend snug against bottom of plot window
-;               /top,pos=[x,y] and /bottom,pos=[x,y] produce same positions.
-;
-;       If LINESTYLE, PSYM, VECTORFONT, THICK, COLORS, or TEXTCOLORS are
-;       supplied as scalars, then the scalar value is set for every line or
-;       symbol in the legend.
-; Outputs:
-;       legend to current plot device
-; OPTIONAL OUTPUT KEYWORDS:
-;       corners = 4-element array, like !p.position, of the normalized
-;         coords for the box (even if box=0): [llx,lly,urx,ury].
-;         Useful for multi-column or multi-line legends, for example,
-;         to make a 2-column legend, you might do the following:
-;           c1_items = ['diamond','asterisk','square']
-;           c1_psym = [4,2,6]
-;           c2_items = ['solid','dashed','dotted']
-;           c2_line = [0,2,1]
-;           legend,c1_items,psym=c1_psym,corners=c1,box=0
-;           legend,c2_items,line=c2_line,corners=c2,box=0,pos=[c1[2],c1[3]]
-;           c = [c1[0]<c2[0],c1[1]<c2[1],c1[2]>c2[2],c1[3]>c2[3]]
-;           plots,[c[0],c[0],c[2],c[2],c[0]],[c[1],c[3],c[3],c[1],c[1]],/norm
-;         Useful also to place the legend.  Here's an automatic way to place
-;         the legend in the lower right corner.  The difficulty is that the
-;         legend's width is unknown until it is plotted.  In this example,
-;         the legend is plotted twice: the first time in the upper left, the
-;         second time in the lower right.
-;           legend,['1','22','333','4444'],linestyle=indgen(4),corners=corners
-;                       ; BOGUS LEGEND---FIRST TIME TO REPORT CORNERS
-;           xydims = [corners[2]-corners[0],corners[3]-corners[1]]
-;                       ; SAVE WIDTH AND HEIGHT
-;           chdim=[!d.x_ch_size/float(!d.x_size),!d.y_ch_size/float(!d.y_size)]
-;                       ; DIMENSIONS OF ONE CHARACTER IN NORMALIZED COORDS
-;           pos = [!x.window[1]-chdim[0]-xydims[0] $
-;                       ,!y.window[0]+chdim[1]+xydims[1]]
-;                       ; CALCULATE POSITION FOR LOWER RIGHT
-;           plot,findgen(10)    ; SIMPLE PLOT; YOU DO WHATEVER YOU WANT HERE.
-;           legend,['1','22','333','4444'],linestyle=indgen(4),pos=pos
-;                       ; REDO THE LEGEND IN LOWER RIGHT CORNER
-;         You can modify the pos calculation to place the legend where you
-;         want.  For example to place it in the upper right:
-;           pos = [!x.window[1]-chdim[0]-xydims[0],!y.window[1]-xydims[1]]
-; Common blocks:
-;       none
-; Procedure:
-;       If keyword help is set, call doc_library to print header.
-;       See notes in the code.  Much of the code deals with placement of the
-;       legend.  The main problem with placement is not being
-;       able to sense the length of a string before it is output.  Some crude
-;       approximations are used for centering.
-; Restrictions:
-;       Here are some things that aren't implemented.
-;       - An orientation keyword would allow lines at angles in the legend.
-;       - An array of usersyms would be nice---simple change.
-;       - An order option to interchange symbols and text might be nice.
-;       - Somebody might like double boxes, e.g., with box = 2.
-;       - Another feature might be a continuous bar with ticks and text.
-;       - There are no guards to avoid writing outside the plot area.
-;       - There is no provision for multi-line text, e.g., '1st line!c2nd line'
-;         Sensing !c would be easy, but !c isn't implemented for PostScript.
-;         A better way might be to simply output the 2nd line as another item
-;         but without any accompanying symbol or linestyle.  A flag to omit
-;         the symbol and linestyle is linestyle[i] = -1.
-;       - There is no ability to make a title line containing any of titles
-;         for the legend, for the symbols, or for the text.
-; Side Effects:
-; Modification history:
-;       write, 24-25 Aug 92, F K Knight (knight@ll.mit.edu)
-;       allow omission of items or omission of both psym and linestyle, add
-;         corners keyword to facilitate multi-column legends, improve place-
-;         ment of symbols and text, add guards for unequal size, 26 Aug 92, FKK
-;       add linestyle(i)=-1 to suppress a single symbol/line, 27 Aug 92, FKK
-;       add keyword vectorfont to allow characters in the sym/line column,
-;         28 Aug 92, FKK
-;       add /top, /bottom, /left, /right keywords for automatic placement at
-;         the four corners of the plot window.  The /right keyword forces
-;         right-to-left printing of menu. 18 Jun 93, FKK
-;       change default position to data coords and add normal, data, and
-;         device keywords, 17 Jan 94, FKK
-;       add /center keyword for positioning, but it is not precise because
-;         text string lengths cannot be known in advance, 17 Jan 94, FKK
-;       add interactive positioning with /position keyword, 17 Jan 94, FKK
-;       allow a legend with just text, no plotting symbols.  This helps in
-;         simply describing a plot or writing assumptions done, 4 Feb 94, FKK
-;       added thick, symsize, and clear keyword Feb 96, W. Landsman HSTX
-;               David Seed, HR Wallingford, d.seed@hrwallingford.co.uk
-;       allow scalar specification of keywords, Mar 96, W. Landsman HSTX
-;       added charthick keyword, June 96, W. Landsman HSTX
-;       Made keyword names  left,right,top,bottom,center longer,
-;                                 Aug 16, 2000, Kim Tolbert
-;       Added ability to have regular text lines in addition to plot legend 
-;       lines in legend.  If linestyle is -99 that item is left-justified.
-;       Previously, only option for no sym/line was linestyle=-1, but then text
-;       was lined up after sym/line column.    10 Oct 2000, Kim Tolbert
-;       Make default value of thick = !P.thick  W. Landsman  Jan. 2001
-;       Don't overwrite existing USERSYM definition  W. Landsman Mar. 2002
-;-
-pro legend, items, BOTTOM_LEGEND=bottom, BOX = box, CENTER_LEGEND=center, $
-    CHARTHICK=charthick, CHARSIZE = charsize, CLEAR = clear, COLORS = colorsi, $
-    CORNERS = corners, DATA=data, DELIMITER=delimiter, DEVICE=device, $
-    FILL=fill, HELP = help, HORIZONTAL=horizontal,LEFT_LEGEND=left, $
-    LINESTYLE=linestylei, MARGIN=margin, NORMAL=normal, NUMBER=number, $
-    POSITION=position,PSPACING=pspacing, PSYM=psymi, RIGHT_LEGEND=right, $
-    SPACING=spacing, SYMSIZE=symsize, TEXTCOLORS=textcolorsi, THICK=thicki, $
-    TOP_LEGEND=top, USERSYM=usersym,  VECTORFONT=vectorfonti, VERTICAL=vertical
-;
-;       =====>> HELP
-;
-on_error,2
-if keyword_set(help) then begin & doc_library,'legend' & return & endif
-;
-;       =====>> SET DEFAULTS FOR SYMBOLS, LINESTYLES, AND ITEMS.
-;
- ni = n_elements(items)
- np = n_elements(psymi)
- nl = n_elements(linestylei)
- nth = n_elements(thicki)
- nv = n_elements(vectorfonti)
- nlpv = max([np,nl,nv])
- n = max([ni,np,nl,nv])                                  ; NUMBER OF ENTRIES
-strn = strtrim(n,2)                                     ; FOR ERROR MESSAGES
-if n eq 0 then message,'No inputs!  For help, type legend,/help.'
-if ni eq 0 then begin
-  items = replicate('',n)                               ; DEFAULT BLANK ARRAY
-endif else begin
-  if size(items,/TNAME) NE 'STRING' then message, $
-      'First parameter must be a string array.  For help, type legend,/help.'
-  if ni ne n then message,'Must have number of items equal to '+strn
-endelse
-symline = (np ne 0) or (nl ne 0)                        ; FLAG TO PLOT SYM/LINE
- if (np ne 0) and (np ne n) and (np NE 1) then message, $
-        'Must have 0, 1 or '+strn+' elements in PSYM array.'
- if (nl ne 0) and (nl ne n) and (nl NE 1) then message, $
-         'Must have 0, 1 or '+strn+' elements in LINESTYLE array.'
- if (nth ne 0) and (nth ne n) and (nth NE 1) then message, $
-         'Must have 0, 1 or '+strn+' elements in THICK array.'
-
- case nl of 
- 0: linestyle = intarr(n)              ;Default = solid
- 1: linestyle = intarr(n)  + linestylei
- else: linestyle = linestylei
- endcase 
- 
- case nth of 
- 0: thick = replicate(!p.thick,n)      ;Default = !P.THICK
- 1: thick = intarr(n) + thicki
- else: thick = thicki
- endcase 
-
- case np of             ;Get symbols
- 0: psym = intarr(n)    ;Default = solid
- 1: psym = intarr(n) + psymi
- else: psym = psymi
- endcase 
-
- case nv of 
- 0: vectorfont = replicate('',n)
- 1: vectorfont = replicate(vectorfonti,n)
- else: vectorfont = vectorfonti
- endcase 
-;
-;       =====>> CHOOSE VERTICAL OR HORIZONTAL ORIENTATION.
-;
-if n_elements(horizontal) eq 0 then begin               ; D=VERTICAL
-  if n_elements(vertical) eq 0 then vertical = 1
-endif else begin
-  if n_elements(vertical) eq 0 then vertical = not horizontal
-endelse
-;
-;       =====>> SET DEFAULTS FOR OTHER OPTIONS.
-;
-if n_elements(box) eq 0 then box = 1
-if n_elements(clear) eq 0 then clear = 0
-
-if n_elements(margin) eq 0 then margin = 0.5
-if n_elements(delimiter) eq 0 then delimiter = ''
-if n_elements(charsize) eq 0 then charsize = !p.charsize
-if n_elements(charthick) eq 0 then charthick = !p.charthick
-if charsize eq 0 then charsize = 1
-if (n_elements (symsize) eq 0) then symsize= charsize + intarr(n)
-if n_elements(number) eq 0 then number = 1
- case N_elements(colorsi) of 
- 0: colors = replicate(!P.color,n)     ;Default is !P.COLOR
- 1: colors = replicate(colorsi,n)
- else: colors = colorsi
- endcase 
-
- case N_elements(textcolorsi) of 
- 0: textcolors = replicate(!P.color,n)      ;Default is !P.COLOR
- 1: textcolors = replicate(textcolorsi,n)
- else: textcolors = textcolorsi
- endcase 
- fill = keyword_set(fill)
-if n_elements(usersym) eq 1 then usersym = 2*[[0,0],[0,1],[1,1],[1,0],[0,0]]-1
-;
-;       =====>> INITIALIZE SPACING
-;
-if n_elements(spacing) eq 0 then spacing = 1.2
-if n_elements(pspacing) eq 0 then pspacing = 3
-xspacing = !d.x_ch_size/float(!d.x_size) * (spacing > charsize)
-yspacing = !d.y_ch_size/float(!d.y_size) * (spacing > charsize)
-ltor = 1                                        ; flag for left-to-right
-if n_elements(left) eq 1 then ltor = left eq 1
-if n_elements(right) eq 1 then ltor = right ne 1
-ttob = 1                                        ; flag for top-to-bottom
-if n_elements(top) eq 1 then ttob = top eq 1
-if n_elements(bottom) eq 1 then ttob = bottom ne 1
-xalign = ltor ne 1                              ; x alignment: 1 or 0
-yalign = -0.5*ttob + 1                          ; y alignment: 0.5 or 1
-xsign = 2*ltor - 1                              ; xspacing direction: 1 or -1
-ysign = 2*ttob - 1                              ; yspacing direction: 1 or -1
-if not ttob then yspacing = -yspacing
-if not ltor then xspacing = -xspacing
-;
-;       =====>> INITIALIZE POSITIONS: FIRST CALCULATE X OFFSET FOR TEXT
-;
-xt = 0
-if nlpv gt 0 then begin                         ; SKIP IF TEXT ITEMS ONLY.
-if vertical then begin                          ; CALC OFFSET FOR TEXT START
-  for i = 0,n-1 do begin
-    if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
-    if psym[i] lt 0 then num = number > 2       ; TO SHOW CONNECTING LINE
-    if psym[i] eq 0 then expand = 1 else expand = 2
-    thisxt = (expand*pspacing*(num-1)*xspacing)
-    if ltor then xt = thisxt > xt else xt = thisxt < xt
-    endfor
-endif   ; NOW xt IS AN X OFFSET TO ALIGN ALL TEXT ENTRIES.
-endif
-;
-;       =====>> INITIALIZE POSITIONS: SECOND LOCATE BORDER
-;
-if !x.window[0] eq !x.window[1] then begin
-  plot,/nodata,xstyle=4,ystyle=4,[0],/noerase
-endif
-;       next line takes care of weirdness with small windows
-pos = [min(!x.window),min(!y.window),max(!x.window),max(!y.window)]
-case n_elements(position) of
- 0: begin
-  if ltor then px = pos[0] else px = pos[2]
-  if ttob then py = pos[3] else py = pos[1]
-  if keyword_set(center) then begin
-    if not keyword_set(right) and not keyword_set(left) then $
-      px = (pos[0] + pos[2])/2. - xt
-    if not keyword_set(top) and not keyword_set(bottom) then $
-      py = (pos[1] + pos[3])/2. + n*yspacing
-    endif
-  position = [px,py] + [xspacing,-yspacing]
-  end
- 1: begin       ; interactive
-  message,/inform,'Place mouse at upper left corner and click any mouse button.'
-  cursor,x,y,/normal
-  position = [x,y]
-  end
- 2: begin       ; convert upper left corner to normal coordinates
-  if keyword_set(data) then $
-    position = convert_coord(position,/to_norm) $
-  else if keyword_set(device) then $
-    position = convert_coord(position,/to_norm,/device) $
-  else if not keyword_set(normal) then $
-    position = convert_coord(position,/to_norm)
-  end
- else: message,'Position keyword can have 0, 1, or 2 elements only. Try legend,/help.'
-endcase
-
-yoff = 0.25*yspacing*ysign                      ; VERT. OFFSET FOR SYM/LINE.
-
-x0 = position[0] + (margin)*xspacing            ; INITIAL X & Y POSITIONS
-y0 = position[1] - margin*yspacing + yalign*yspacing    ; WELL, THIS WORKS!
-;
-;       =====>> OUTPUT TEXT FOR LEGEND, ITEM BY ITEM.
-;       =====>> FOR EACH ITEM, PLACE SYM/LINE, THEN DELIMITER,
-;       =====>> THEN TEXT---UPDATING X & Y POSITIONS EACH TIME.
-;       =====>> THERE ARE A NUMBER OF EXCEPTIONS DONE WITH IF STATEMENTS.
-;
-for iclr = 0,clear do begin
-  y = y0                                                ; STARTING X & Y POSITIONS
-  x = x0
-  if ltor then xend = 0 else xend = 1           ; SAVED WIDTH FOR DRAWING BOX
-
- if ttob then ii = [0,n-1,1] else ii = [n-1,0,-1]
- for i = ii[0],ii[1],ii[2] do begin
-  if vertical then x = x0 else y = y0           ; RESET EITHER X OR Y
-  x = x + xspacing                              ; UPDATE X & Y POSITIONS
-  y = y - yspacing
-  if nlpv eq 0 then goto,TEXT_ONLY              ; FLAG FOR TEXT ONLY
-  if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
-  if psym[i] lt 0 then num = number > 2         ; TO SHOW CONNECTING LINE
-  if psym[i] eq 0 then expand = 1 else expand = 2
-  xp = x + expand*pspacing*indgen(num)*xspacing
-  if (psym[i] gt 0) and (num eq 1) and vertical then xp = x + xt/2.
-  yp = y + intarr(num)
-  if vectorfont[i] eq '' then yp = yp + yoff
-  if psym[i] eq 0 then begin
-    xp = [min(xp),max(xp)]                      ; TO EXPOSE LINESTYLES
-    yp = [min(yp),max(yp)]                      ; DITTO
-    endif
-  if (psym[i] eq 8) and (N_elements(usersym) GT 1) then $
-                usersym,usersym,fill=fill,color=colors[i]
-;; extra by djseed .. psym=88 means use the already defined usersymbol
- if psym[i] eq 88 then psym[i] =8
-  if vectorfont[i] ne '' then begin
-;    if (num eq 1) and vertical then xp = x + xt/2      ; IF 1, CENTERED.
-    xyouts,xp,yp,vectorfont[i],width=width,color=colors[i] $
-      ,size=charsize,align=xalign,charthick = charthick,/norm
-    xt = xt > width
-    xp = xp + width/2.
-  endif else begin
-    if symline and (linestyle[i] ge 0) then plots,xp,yp,color=colors[i] $
-      ,/normal,linestyle=linestyle[i],psym=psym[i],symsize=symsize[i], $
-      thick=thick[i]
-  endelse
-
-  if vertical then x = x + xt else if ltor then x = max(xp) else x = min(xp)
-  if symline then x = x + xspacing
-  TEXT_ONLY:
-  if vertical and (vectorfont[i] eq '') and symline and (linestyle[i] eq -99) then x=x0 + xspacing
-  xyouts,x,y,delimiter,width=width,/norm,color=textcolors[i], $
-         size=charsize,align=xalign,charthick = charthick
-  x = x + width*xsign
-  if width ne 0 then x = x + 0.5*xspacing
-  xyouts,x,y,items[i],width=width,/norm,color=textcolors[i],size=charsize, $
-             align=xalign,charthick=charthick
-  x = x + width*xsign
-  if not vertical and (i lt (n-1)) then x = x+2*xspacing; ADD INTER-ITEM SPACE
-  xfinal = (x + xspacing*margin)
-  if ltor then xend = xfinal > xend else xend = xfinal < xend   ; UPDATE END X
- endfor
-
- if (iclr lt clear ) then begin
-;       =====>> CLEAR AREA
-        x = position[0]
-        y = position[1]
-        if vertical then bottom = n else bottom = 1
-        ywidth = - (2*margin+bottom-0.5)*yspacing
-        corners = [x,y+ywidth,xend,y]
-        polyfill,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm,color=-1
-;       plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],thick=2
- endif else begin
-
-;
-;       =====>> OUTPUT BORDER
-;
-        x = position[0]
-        y = position[1]
-        if vertical then bottom = n else bottom = 1
-        ywidth = - (2*margin+bottom-0.5)*yspacing
-        corners = [x,y+ywidth,xend,y]
-        if box then plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm
-        return
- endelse
-endfor
-
-end
-
diff --git a/examples/IsothermalPotential/GravityOnly/run.sh b/examples/IsothermalPotential/GravityOnly/run.sh
index 329d04d30a6ccba7036f3802ff184c854ed00761..f6adfcececf4923485c0deabd97e9af9a6f64b05 100755
--- a/examples/IsothermalPotential/GravityOnly/run.sh
+++ b/examples/IsothermalPotential/GravityOnly/run.sh
@@ -7,4 +7,4 @@ then
     python makeIC.py 1000 1
 fi
 
-../../swift -g -t 2 isothermal.yml
+../../swift -g -t 2 isothermal.yml 2>&1 | tee output.log
diff --git a/examples/KelvinHelmoltz_2D/kelvinHelmholtz.yml b/examples/KelvinHelmholtz_2D/kelvinHelmholtz.yml
similarity index 100%
rename from examples/KelvinHelmoltz_2D/kelvinHelmholtz.yml
rename to examples/KelvinHelmholtz_2D/kelvinHelmholtz.yml
diff --git a/examples/KelvinHelmoltz_2D/makeIC.py b/examples/KelvinHelmholtz_2D/makeIC.py
similarity index 100%
rename from examples/KelvinHelmoltz_2D/makeIC.py
rename to examples/KelvinHelmholtz_2D/makeIC.py
diff --git a/examples/KelvinHelmoltz_2D/plotSolution.py b/examples/KelvinHelmholtz_2D/plotSolution.py
similarity index 100%
rename from examples/KelvinHelmoltz_2D/plotSolution.py
rename to examples/KelvinHelmholtz_2D/plotSolution.py
diff --git a/examples/KelvinHelmoltz_2D/run.sh b/examples/KelvinHelmholtz_2D/run.sh
similarity index 82%
rename from examples/KelvinHelmoltz_2D/run.sh
rename to examples/KelvinHelmholtz_2D/run.sh
index 4899ca89bc7bbbf72d15d6ecd3961c146a9c9821..3d83bcc6dd7c226885ed83c3188f3177c4807154 100755
--- a/examples/KelvinHelmoltz_2D/run.sh
+++ b/examples/KelvinHelmholtz_2D/run.sh
@@ -8,7 +8,7 @@ then
 fi
 
 # Run SWIFT
-../swift -s -t 1 kelvinHelmholtz.yml
+../swift -s -t 1 kelvinHelmholtz.yml 2>&1 | tee output.log
 
 # Plot the solution
 python plotSolution.py 6
diff --git a/examples/MultiTypes/run.sh b/examples/MultiTypes/run.sh
index 9782a27fe18a7fac94cc2a12353f81abb9efb06b..57465ce0ba6dde3988359df990f2a93323dbc617 100755
--- a/examples/MultiTypes/run.sh
+++ b/examples/MultiTypes/run.sh
@@ -7,4 +7,4 @@ then
     python makeIC.py 50 60
 fi
 
-../swift -s -g -t 16 multiTypes.yml
+../swift -s -g -t 16 multiTypes.yml 2>&1 | tee output.log
diff --git a/examples/PerturbedBox_3D/run.sh b/examples/PerturbedBox_3D/run.sh
index c15011ca526efb7250e93dc8cb002e49647b82b6..e20bff52d18322ce377fb589900fd9e13eefe64d 100755
--- a/examples/PerturbedBox_3D/run.sh
+++ b/examples/PerturbedBox_3D/run.sh
@@ -7,4 +7,4 @@ then
     python makeIC.py 50
 fi
 
-../swift -s -t 16 perturbedBox.yml
+../swift -s -t 16 perturbedBox.yml 2>&1 | tee output.log
diff --git a/examples/SedovBlast_1D/run.sh b/examples/SedovBlast_1D/run.sh
index 911b595a604cae2db72c6cab070189848e677010..4b9a84f069673bd6def3b96faec71b9d4fdd0dda 100755
--- a/examples/SedovBlast_1D/run.sh
+++ b/examples/SedovBlast_1D/run.sh
@@ -8,7 +8,7 @@ then
 fi
 
 # Run SWIFT
-../swift -s -t 1 sedov.yml
+../swift -s -t 1 sedov.yml 2>&1 | tee output.log
 
 # Plot the solution
 python plotSolution.py 5
diff --git a/examples/SedovBlast_2D/run.sh b/examples/SedovBlast_2D/run.sh
index b6f6f73db3ff281492d2abc21eb5ded67e528516..a32c8f0d6f3116d5486fe1bd086bf8df49d06020 100755
--- a/examples/SedovBlast_2D/run.sh
+++ b/examples/SedovBlast_2D/run.sh
@@ -13,7 +13,7 @@ then
 fi
 
 # Run SWIFT
-../swift -s -t 1 sedov.yml
+../swift -s -t 1 sedov.yml 2>&1 | tee output.log
 
 # Plot the solution
 python plotSolution.py 5
diff --git a/examples/SedovBlast_3D/run.sh b/examples/SedovBlast_3D/run.sh
index afc9ff257eff334d700dbd139c7ffbf0225fbf7b..00d5e5b91c31e64f824a3d2a28c8e1a126684a74 100755
--- a/examples/SedovBlast_3D/run.sh
+++ b/examples/SedovBlast_3D/run.sh
@@ -13,7 +13,7 @@ then
 fi
 
 # Run SWIFT
-../swift -s -t 4 sedov.yml
+../swift -s -t 4 sedov.yml 2>&1 | tee output.log
 
 # Plot the solution
 python plotSolution.py 5
diff --git a/examples/SodShock_1D/run.sh b/examples/SodShock_1D/run.sh
index e3ac218c56caa81d0e7a6817e03a8db20bb575d5..4be4254baa4a87b105a5f3c1bfbf9059348a1e9e 100755
--- a/examples/SodShock_1D/run.sh
+++ b/examples/SodShock_1D/run.sh
@@ -8,7 +8,7 @@ then
 fi
 
 # Run SWIFT
-../swift -s -t 1 sodShock.yml
+../swift -s -t 1 sodShock.yml 2>&1 | tee output.log
 
 # Plot the result
-python plotSolution.py 1
+python plotSolution.py 1 
diff --git a/examples/SodShock_2D/run.sh b/examples/SodShock_2D/run.sh
index 3d1f4beb3ac4c4becfda24ed39390efec8074da7..9e6bbfdf1c0a7c206ce6966fdca7b20a28047dd8 100755
--- a/examples/SodShock_2D/run.sh
+++ b/examples/SodShock_2D/run.sh
@@ -13,6 +13,6 @@ then
 fi
 
 # Run SWIFT
-../swift -s -t 1 sodShock.yml
+../swift -s -t 1 sodShock.yml 2>&1 | tee output.log
 
 python plotSolution.py 1
diff --git a/examples/SodShock_3D/run.sh b/examples/SodShock_3D/run.sh
index 0f9e63be334475d98196189c49b95fc46982704a..8ed85baf73425b75f402c491a3c66785f6c6fce0 100755
--- a/examples/SodShock_3D/run.sh
+++ b/examples/SodShock_3D/run.sh
@@ -13,6 +13,6 @@ then
 fi
 
 # Run SWIFT
-../swift -s -t 4 sodShock.yml
+../swift -s -t 4 sodShock.yml 2>&1 | tee output.log
 
 python plotSolution.py 1
diff --git a/examples/SquareTest_2D/run.sh b/examples/SquareTest_2D/run.sh
index 21f5ee9093ca2b50830c46bb2117d8138d23ee94..242f1b8b729979b026cfe31002e84cd9ef741129 100755
--- a/examples/SquareTest_2D/run.sh
+++ b/examples/SquareTest_2D/run.sh
@@ -8,7 +8,7 @@ then
 fi
 
 # Run SWIFT
-../swift -s -t 1 square.yml
+../swift -s -t 1 square.yml 2>&1 | tee output.log
 
 # Plot the solution
 python plotSolution.py 40
diff --git a/examples/UniformBox_2D/run.sh b/examples/UniformBox_2D/run.sh
index fbabedb1cd564e0831228591ce32d9a73fe08d9a..ee3ef109968a65e2437ea17b42013266195d3314 100755
--- a/examples/UniformBox_2D/run.sh
+++ b/examples/UniformBox_2D/run.sh
@@ -7,4 +7,4 @@ then
     python makeIC.py 100
 fi
 
-../swift -s -t 16 uniformPlane.yml
+../swift -s -t 16 uniformPlane.yml 2>&1 | tee output.log
diff --git a/examples/UniformBox_3D/run.sh b/examples/UniformBox_3D/run.sh
index 0cb0a505915be47bafb99ed7531685bfeb3dc829..08891cdd08fccf8f43089951e94dddb33e162030 100755
--- a/examples/UniformBox_3D/run.sh
+++ b/examples/UniformBox_3D/run.sh
@@ -7,4 +7,4 @@ then
     python makeIC.py 100
 fi
 
-../swift -s -t 16 uniformBox.yml
+../swift -s -t 16 uniformBox.yml 2>&1 | tee output.log
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index b89b4f70aea1c29e2df4df3615e34552e4ea3f8d..96543c933048ed9404d5602f1a1986b7cc34ed28 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -107,9 +107,10 @@ void writeSPHflavour(hid_t h_grpsph) {
 
   /* Viscosity and thermal conduction */
   writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
-                   "(No treatment) Legacy Gadget-2 as in Springel (2005)");
-  writeAttribute_s(h_grpsph, "Viscosity Model",
-                   "Legacy Gadget-2 as in Springel (2005)");
+                   "(No treatment) as in Springel (2005)");
+  writeAttribute_s(
+      h_grpsph, "Viscosity Model",
+      "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
   writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
   writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
 }
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 842a028b9031b3c1cf36e4735089d05e038abdca..a5d73aad02372aa22b840c2cb0d3100cd439e75d 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -16,10 +16,29 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
+#ifndef SWIFT_MINIMAL_HYDRO_H
+#define SWIFT_MINIMAL_HYDRO_H
+
+/**
+ * @file Minimal/hydro.h
+ * @brief Minimal conservative implementation of SPH (Non-neighbour loop
+ * equations)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
+ */
 
 #include "adiabatic_index.h"
-#include "approx_math.h"
+#include "dimension.h"
 #include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "kernel_hydro.h"
 
 /**
  * @brief Returns the internal energy of a particle
@@ -34,7 +53,7 @@
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
     const struct part *restrict p, float dt) {
 
-  return p->u;
+  return p->u + p->u_dt * dt;
 }
 
 /**
@@ -46,7 +65,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
 __attribute__((always_inline)) INLINE static float hydro_get_pressure(
     const struct part *restrict p, float dt) {
 
-  return gas_pressure_from_internal_energy(p->rho, p->u);
+  const float u = p->u + p->u_dt * dt;
+
+  return gas_pressure_from_internal_energy(p->rho, u);
 }
 
 /**
@@ -62,7 +83,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
 __attribute__((always_inline)) INLINE static float hydro_get_entropy(
     const struct part *restrict p, float dt) {
 
-  return gas_entropy_from_internal_energy(p->rho, p->u);
+  const float u = p->u + p->u_dt * dt;
+
+  return gas_entropy_from_internal_energy(p->rho, u);
 }
 
 /**
@@ -74,7 +97,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
 __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
     const struct part *restrict p, float dt) {
 
-  return gas_soundspeed_from_internal_energy(p->rho, p->u);
+  const float u = p->u + p->u_dt * dt;
+
+  return gas_soundspeed_from_internal_energy(p->rho, u);
 }
 
 /**
@@ -150,7 +175,6 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
   xp->v_full[0] = p->v[0];
   xp->v_full[1] = p->v[1];
   xp->v_full[2] = p->v[2];
-  xp->u_full = p->u;
 }
 
 /**
@@ -164,6 +188,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
  */
 __attribute__((always_inline)) INLINE static void hydro_init_part(
     struct part *restrict p) {
+
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
@@ -227,7 +252,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp, int ti_current,
     double timeBase) {
 
-  p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
+  const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
+  const float pressure = hydro_get_pressure(p, half_dt);
+
+  p->force.pressure = pressure;
 }
 
 /**
@@ -268,10 +296,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     struct part *restrict p, const struct xpart *restrict xp, int t0, int t1,
     double timeBase) {
 
-  p->u = xp->u_full;
-
-  /* Need to recompute the pressure as well */
-  p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
+  /* Drift the pressure */
+  const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
+  p->force.pressure = hydro_get_pressure(p, dt_entr);
 }
 
 /**
@@ -304,11 +331,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     struct part *restrict p, struct xpart *restrict xp, float dt,
     float half_dt) {
 
-  /* Kick in momentum space */
-  xp->u_full += p->u_dt * dt;
+  /* Do not decrease the energy by more than a factor of 2*/
+  const float u_change = p->u_dt * dt;
+  if (u_change > -0.5f * p->u)
+    p->u += u_change;
+  else
+    p->u *= 0.5f;
 
-  /* Get the predicted internal energy */
-  p->u = xp->u_full - half_dt * p->u_dt;
+  /* Do not 'overcool' when timestep increases */
+  if (p->u + p->u_dt * half_dt < 0.5f * p->u) p->u_dt = -0.5f * p->u / half_dt;
 }
 
 /**
@@ -323,3 +354,5 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
  */
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
     struct part *restrict p) {}
+
+#endif /* SWIFT_MINIMAL_HYDRO_H */
diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index 7ab348d255173a33fc29e482eb6032c03c82d9eb..0a6e50da6051fd961f4bf35f32caee311888a13e 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -16,18 +16,36 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
+#ifndef SWIFT_MINIMAL_HYDRO_DEBUG_H
+#define SWIFT_MINIMAL_HYDRO_DEBUG_H
+
+/**
+ * @file Minimal/hydro_debug.h
+ * @brief Minimal conservative implementation of SPH (Debugging routines)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
+ */
 
 __attribute__((always_inline)) INLINE static void hydro_debug_particle(
     const struct part* p, const struct xpart* xp) {
   printf(
       "x=[%.3e,%.3e,%.3e], "
       "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
-      "u_full=%.3e, u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
-      "h=%.3e, dh/dt=%.3e "
-      "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%d, t_end=%d\n",
+      "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
+      "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, "
+      "t_begin=%d, t_end=%d\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
-      xp->u_full, p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h,
-      p->force.h_dt, (int)p->density.wcount, p->mass, p->rho_dh, p->rho,
-      p->ti_begin, p->ti_end);
+      p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->force.h_dt,
+      (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
+      p->ti_end);
 }
+
+#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 9eda45760248ebc0dde18e503abe0b1288659518..a8a855d9db81f6927c1d8b45410a57d50a8366de 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -16,18 +16,25 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef SWIFT_RUNNER_IACT_MINIMAL_H
-#define SWIFT_RUNNER_IACT_MINIMAL_H
-
-#include "adiabatic_index.h"
+#ifndef SWIFT_MINIMAL_HYDRO_IACT_H
+#define SWIFT_MINIMAL_HYDRO_IACT_H
 
 /**
- * @brief Minimal conservative implementation of SPH
+ * @file Minimal/hydro_iact.h
+ * @brief Minimal conservative implementation of SPH (Neighbour loop equations)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
  *
- * The thermal variable is the internal energy (u). No viscosity nor
- * thermal conduction terms are implemented.
+ * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
  */
 
+#include "adiabatic_index.h"
+
 /**
  * @brief Density loop
  */
@@ -115,6 +122,8 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
 __attribute__((always_inline)) INLINE static void runner_iact_force(
     float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
 
+  const float fac_mu = 1.f; /* Will change with cosmological integration */
+
   const float r = sqrtf(r2);
   const float r_inv = 1.0f / r;
 
@@ -143,34 +152,60 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float wj_dr = hjd_inv * wj_dx;
 
   /* Compute gradient terms */
-  const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
-  const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
+  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
 
   /* Compute dv dot r. */
   const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
                      (pi->v[1] - pj->v[1]) * dx[1] +
                      (pi->v[2] - pj->v[2]) * dx[2];
 
-  /* Compute sound speeds */
+  /* Are the particles moving towards each others ? */
+  const float omega_ij = fminf(dvdr, 0.f);
+  const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+  /* Compute sound speeds and signal velocity */
   const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
   const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
-  const float v_sig = ci + cj;
+  const float v_sig = ci + cj - 3.f * mu_ij;
+
+  /* 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;
+
+  /* Convolve with the kernel */
+  const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
 
   /* SPH acceleration term */
-  const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
+  const float sph_acc_term =
+      (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+
+  /* Assemble the acceleration */
+  const float acc = sph_acc_term + visc_acc_term;
 
   /* Use the force Luke ! */
-  pi->a_hydro[0] -= mj * sph_term * dx[0];
-  pi->a_hydro[1] -= mj * sph_term * dx[1];
-  pi->a_hydro[2] -= mj * sph_term * dx[2];
+  pi->a_hydro[0] -= mj * acc * dx[0];
+  pi->a_hydro[1] -= mj * acc * dx[1];
+  pi->a_hydro[2] -= mj * acc * dx[2];
 
-  pj->a_hydro[0] += mi * sph_term * dx[0];
-  pj->a_hydro[1] += mi * sph_term * dx[1];
-  pj->a_hydro[2] += mi * sph_term * dx[2];
+  pj->a_hydro[0] += mi * acc * dx[0];
+  pj->a_hydro[1] += mi * acc * dx[1];
+  pj->a_hydro[2] += mi * acc * dx[2];
 
   /* Get the time derivative for u. */
-  pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
-  pj->u_dt += P_over_rho_j * mi * dvdr * r_inv * wj_dr;
+  const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
+  const float sph_du_term_j = P_over_rho2_j * dvdr * r_inv * wj_dr;
+
+  /* Viscosity term */
+  const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+
+  /* Assemble the energy equation term */
+  const float du_dt_i = sph_du_term_i + visc_du_term;
+  const float du_dt_j = sph_du_term_j + visc_du_term;
+
+  /* Internal energy time derivatibe */
+  pi->u_dt += du_dt_i * mj;
+  pj->u_dt += du_dt_j * mi;
 
   /* Get the time derivative for h. */
   pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
@@ -198,6 +233,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
 __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
     float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
 
+  const float fac_mu = 1.f; /* Will change with cosmological integration */
+
   const float r = sqrtf(r2);
   const float r_inv = 1.0f / r;
 
@@ -226,29 +263,53 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float wj_dr = hjd_inv * wj_dx;
 
   /* Compute gradient terms */
-  const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
-  const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
+  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
 
   /* Compute dv dot r. */
   const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
                      (pi->v[1] - pj->v[1]) * dx[1] +
                      (pi->v[2] - pj->v[2]) * dx[2];
 
-  /* Compute sound speeds */
+  /* Are the particles moving towards each others ? */
+  const float omega_ij = fminf(dvdr, 0.f);
+  const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+  /* Compute sound speeds and signal velocity */
   const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
   const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
-  const float v_sig = ci + cj;
+  const float v_sig = ci + cj - 3.f * mu_ij;
+
+  /* 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;
+
+  /* Convolve with the kernel */
+  const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
 
   /* SPH acceleration term */
-  const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
+  const float sph_acc_term =
+      (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+
+  /* Assemble the acceleration */
+  const float acc = sph_acc_term + visc_acc_term;
 
   /* Use the force Luke ! */
-  pi->a_hydro[0] -= mj * sph_term * dx[0];
-  pi->a_hydro[1] -= mj * sph_term * dx[1];
-  pi->a_hydro[2] -= mj * sph_term * dx[2];
+  pi->a_hydro[0] -= mj * acc * dx[0];
+  pi->a_hydro[1] -= mj * acc * dx[1];
+  pi->a_hydro[2] -= mj * acc * dx[2];
 
   /* Get the time derivative for u. */
-  pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
+  const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
+
+  /* Viscosity term */
+  const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+
+  /* Assemble the energy equation term */
+  const float du_dt_i = sph_du_term_i + visc_du_term;
+
+  /* Internal energy time derivatibe */
+  pi->u_dt += du_dt_i * mj;
 
   /* Get the time derivative for h. */
   pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
@@ -268,4 +329,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
       "function does not exist yet!");
 }
 
-#endif /* SWIFT_RUNNER_IACT_MINIMAL_H */
+#endif /* SWIFT_MINIMAL_HYDRO_IACT_H */
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 1a5b2938354fb8ee3ef936f86ff442c731e0d4c7..01a75b17fd5577cfcfb48d3afac22579f30fcf7a 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -16,7 +16,25 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
+#ifndef SWIFT_MINIMAL_HYDRO_IO_H
+#define SWIFT_MINIMAL_HYDRO_IO_H
 
+/**
+ * @file Minimal/hydro_io.h
+ * @brief Minimal conservative implementation of SPH (i/o routines)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
+ */
+
+#include "adiabatic_index.h"
+#include "hydro.h"
 #include "io_properties.h"
 #include "kernel_hydro.h"
 
@@ -51,6 +69,16 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
                                 UNIT_CONV_DENSITY, parts, rho);
 }
 
+float convert_S(struct engine* e, struct part* p) {
+
+  return hydro_get_entropy(p, 0);
+}
+
+float convert_P(struct engine* e, struct part* p) {
+
+  return hydro_get_pressure(p, 0);
+}
+
 /**
  * @brief Specifies which particle fields to write to a dataset
  *
@@ -61,7 +89,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
 void hydro_write_particles(struct part* parts, struct io_props* list,
                            int* num_fields) {
 
-  *num_fields = 8;
+  *num_fields = 10;
 
   /* List what we want to write */
   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -80,6 +108,11 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
                                  UNIT_CONV_ACCELERATION, parts, a_hydro);
   list[7] =
       io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+  list[8] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+                                              UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+                                              parts, rho, convert_S);
+  list[9] = io_make_output_field_convert_part(
+      "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, rho, convert_P);
 }
 
 /**
@@ -90,8 +123,9 @@ void writeSPHflavour(hid_t h_grpsph) {
 
   /* Viscosity and thermal conduction */
   /* Nothing in this minimal model... */
-  writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No model");
-  writeAttribute_s(h_grpsph, "Viscosity Model", "No model");
+  writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
+  writeAttribute_s(h_grpsph, "Viscosity Model",
+                   "Minimal treatment as in Monaghan (1992)");
 
   /* Time integration properties */
   writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
@@ -104,3 +138,5 @@ void writeSPHflavour(hid_t h_grpsph) {
  * @return 1 if entropy is in 'internal energy', 0 otherwise.
  */
 int writeEntropyFlag() { return 0; }
+
+#endif /* SWIFT_MINIMAL_HYDRO_IO_H */
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index f4feb2311ad2bfd8c7227f3b072c5971724a5815..ad65f8b44fc67f4aae6470246cbab91bc3710007 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -16,6 +16,22 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
+#ifndef SWIFT_MINIMAL_HYDRO_PART_H
+#define SWIFT_MINIMAL_HYDRO_PART_H
+
+/**
+ * @file Minimal/hydro_part.h
+ * @brief Minimal conservative implementation of SPH (Particle definition)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
+ */
 
 /**
  * @brief Particle fields not needed during the SPH loops over neighbours.
@@ -31,8 +47,6 @@ struct xpart {
 
   float v_full[3]; /*!< Velocity at the last full step. */
 
-  float u_full; /*!< Thermal energy at the last full step. */
-
 } __attribute__((aligned(xpart_align)));
 
 /**
@@ -107,3 +121,5 @@ struct part {
   struct gpart* gpart; /*!< Pointer to corresponding gravity part. */
 
 } __attribute__((aligned(part_align)));
+
+#endif /* SWIFT_MINIMAL_HYDRO_PART_H */
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index c425ecce6de38b32645b367886e039bb950df68d..0e4eaf150c764e27156747bb01db20a71f03c7b5 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -70,7 +70,7 @@ void hydro_props_print(const struct hydro_props *p) {
   message(
       "Hydrodynamic integration: Max change of volume: %.2f "
       "(max|dlog(h)/dt|=%f).",
-      powf(expf(p->log_max_h_change), hydro_dimension), p->log_max_h_change);
+      pow_dimension(expf(p->log_max_h_change)), p->log_max_h_change);
 
   if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
     message("Maximal iterations in ghost task set to %d (default is %d)",
@@ -90,8 +90,8 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
   writeAttribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
   writeAttribute_f(h_grpsph, "Volume log(max(delta h))", p->log_max_h_change);
   writeAttribute_f(h_grpsph, "Volume max change time-step",
-                   powf(expf(p->log_max_h_change), 3.f));
-  writeAttribute_f(h_grpsph, "Max ghost iterations",
+                   pow_dimension(expf(p->log_max_h_change)));
+  writeAttribute_i(h_grpsph, "Max ghost iterations",
                    p->max_smoothing_iterations);
 }
 #endif