Merge branch 'master' into doself2-vectorisation

Conflicts: tests/Makefile.am tests/test125cellsPerturbed.sh.in tests/tolerance_125_perturbed.dat

Merge branch 'master' into doself2-vectorisation
Conflicts: tests/Makefile.am tests/test125cellsPerturbed.sh.in tests/tolerance_125_perturbed.dat
9ebfdd79 · James Willis · 5c9fd282 · f8bb52df · 9ebfdd79 · 9ebfdd79
Commit 9ebfdd79 authored May 21, 2017 by James Willis
--- a/.gitignore
+++ b/.gitignore
@@ -47,6 +47,8 @@ tests/brute_force_27_perturbed.dat
 tests/swift_dopair_27_perturbed.dat
 tests/brute_force_125_standard.dat
 tests/swift_dopair_125_standard.dat
+tests/brute_force_125_perturbed.dat
+tests/swift_dopair_125_perturbed.dat
 tests/testGreetings
 tests/testReading
 tests/input.hdf5
@@ -65,6 +67,7 @@ tests/parser_output.yml
 tests/test27cells.sh
 tests/test27cellsPerturbed.sh
 tests/test125cells.sh
+tests/test125cellsPerturbed.sh
 tests/testPair.sh
 tests/testPairPerturbed.sh
 tests/testParser.sh

--- a/configure.ac
+++ b/configure.ac
@@ -853,6 +853,7 @@ AC_CONFIG_FILES([tests/testPairPerturbed.sh], [chmod +x tests/testPairPerturbed.
 AC_CONFIG_FILES([tests/test27cells.sh], [chmod +x tests/test27cells.sh])
 AC_CONFIG_FILES([tests/test27cellsPerturbed.sh], [chmod +x tests/test27cellsPerturbed.sh])
 AC_CONFIG_FILES([tests/test125cells.sh], [chmod +x tests/test125cells.sh])
+AC_CONFIG_FILES([tests/test125cellsPerturbed.sh], [chmod +x tests/test125cellsPerturbed.sh])
 AC_CONFIG_FILES([tests/testParser.sh], [chmod +x tests/testParser.sh])

 # Save the compilation options

--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -31,8 +31,6 @@ Gravity:
  eta:                   0.025    # Constant dimensionless multiplier for time integration.
  theta:                 0.7      # Opening angle (Multipole acceptance criterion)
  epsilon:               0.0001   # Softening length (in internal units).
-  a_smooth:              1000.
-  r_cut:                 4.
  
 # Parameters for the hydrodynamics scheme
 SPH:

--- a/examples/HydrostaticHalo/density_profile.py
+++ b/examples/HydrostaticHalo/density_profile.py
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
 import numpy as np
 import h5py as h5
-import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
 import sys

 #for the plotting
@@ -46,7 +67,8 @@ for i in range(n_snaps):
    f = h5.File(filename,'r')
    coords_dset = f["PartType0/Coordinates"]
    coords = np.array(coords_dset)
-#translate coords by centre of box
+
+    #translate coords by centre of box
    header = f["Header"]
    snap_time = header.attrs["Time"]
    snap_time_cgs = snap_time * unit_time_cgs
@@ -63,58 +85,46 @@ for i in range(n_snaps):
    bin_width = bin_edges[1] - bin_edges[0]
    hist = np.histogram(r,bins = bin_edges)[0] # number of particles in each bin

-#find the mass in each radial bin

+    #find the mass in each radial bin
    mass_dset = f["PartType0/Masses"]
-#mass of each particles should be equal
+
+    #mass of each particles should be equal
    part_mass = np.array(mass_dset)[0]
    part_mass_cgs = part_mass * unit_mass_cgs
    part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 

    mass_hist = hist * part_mass_over_virial_mass
    radial_bin_mids = np.linspace(bin_width/2.,max_r - bin_width/2.,n_radial_bins)
-#volume in each radial bin
+
+    #volume in each radial bin
    volume = 4.*np.pi * radial_bin_mids**2 * bin_width

-#now divide hist by the volume so we have a density in each bin

+    #now divide hist by the volume so we have a density in each bin
    density = mass_hist / volume

-    ##read the densities
-
-    # density_dset = f["PartType0/Density"]
-    # density = np.array(density_dset)
-    # density_cgs = density * unit_mass_cgs / unit_length_cgs**3
-    # rho = density_cgs * r_vir_cgs**3 / M_vir_cgs
-
    t = np.linspace(10./n_radial_bins,10.0,1000)
    rho_analytic = t**(-2)/(4.*np.pi)

-    #calculate cooling radius
-
-    #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)

-    #initial analytic density profile
-    
+    #initial analytic density profile    
    if (i == 0):
        r_0 = radial_bin_mids[0]
        rho_0 = density[0]
-
        rho_analytic_init = rho_0 * (radial_bin_mids/r_0)**(-2)
-    plt.plot(radial_bin_mids,density/rho_analytic_init,'ko',label = "Average density of shell")
-    #plt.plot(t,rho_analytic,label = "Initial analytic density profile"
-    plt.xlabel(r"$r / r_{vir}$")
-    plt.ylabel(r"$\rho / \rho_{init})$")
-    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    #plt.ylim((1.e-2,1.e1))
-    #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,20),'r',label = "Cooling radius")
-    plt.xlim((radial_bin_mids[0],max_r))
-    plt.ylim((0,20))
-    plt.plot((0,max_r),(1,1))
-    #plt.xscale('log')
-    #plt.yscale('log')
-    plt.legend(loc = "upper right")
+
+    figure()
+    plot(radial_bin_mids,density/rho_analytic_init,'ko',label = "Average density of shell")
+    #plot(t,rho_analytic,label = "Initial analytic density profile")
+    xlabel(r"$r / r_{vir}$")
+    ylabel(r"$\rho / \rho_{init}$")
+    title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
+    xlim((radial_bin_mids[0],max_r))
+    ylim((0,2))
+    plot((0,max_r),(1,1))
+    legend(loc = "upper right")
    plot_filename = "./plots/density_profile/density_profile_%03d.png" %i
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
+    savefig(plot_filename,format = "png")
+    close()

--- a/examples/HydrostaticHalo/internal_energy_profile.py
+++ b/examples/HydrostaticHalo/internal_energy_profile.py
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
 import numpy as np
 import h5py as h5
-import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
 import sys

 def do_binning(x,y,x_bin_edges):
@@ -48,8 +69,6 @@ unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"]
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
 v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
-#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
-#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
 box_centre = np.array(header.attrs["BoxSize"])
@@ -64,7 +83,8 @@ for i in range(n_snaps):
    f = h5.File(filename,'r')
    coords_dset = f["PartType0/Coordinates"]
    coords = np.array(coords_dset)
-#translate coords by centre of box
+
+    #translate coords by centre of box
    header = f["Header"]
    snap_time = header.attrs["Time"]
    snap_time_cgs = snap_time * unit_time_cgs
@@ -75,11 +95,11 @@ for i in range(n_snaps):
    radius_cgs = radius*unit_length_cgs
    radius_over_virial_radius = radius_cgs / r_vir_cgs

-#get the internal energies
+    #get the internal energies
    u_dset = f["PartType0/InternalEnergy"]
    u = np.array(u_dset)

-#make dimensionless
+    #make dimensionless
    u /= v_c**2/(2. * (gamma - 1.))
    r = radius_over_virial_radius

@@ -90,21 +110,16 @@ for i in range(n_snaps):
    radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins) 
    binned_u = u_totals / hist

-    #calculate cooling radius
-
-    #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
-
-    plt.plot(radial_bin_mids,binned_u,'ko',label = "Numerical solution")
-    #plt.plot((0,1),(1,1),label = "Analytic Solution")
-    #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,2),'r',label = "Cooling radius")
-    plt.legend(loc = "lower right")
-    plt.xlabel(r"$r / r_{vir}$")
-    plt.ylabel(r"$u / (v_c^2 / (2(\gamma - 1)) $")
-    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    plt.ylim((0,2))
+    figure()
+    plot(radial_bin_mids,binned_u,'ko',label = "Numerical solution")
+    legend(loc = "lower right")
+    xlabel(r"$r / r_{vir}$")
+    ylabel(r"$u / (v_c^2 / (2(\gamma - 1)) $")
+    title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
+    ylim((0,2))
    plot_filename = "./plots/internal_energy/internal_energy_profile_%03d.png" %i
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
+    savefig(plot_filename,format = "png")
+    close()


        

--- a/examples/HydrostaticHalo/run.sh
+++ b/examples/HydrostaticHalo/run.sh
 #!/bin/bash

 # Generate the initial conditions if they are not present.
-echo "Generating initial conditions for the isothermal potential box example..."
-python makeIC.py 100000
+if [ ! -e Hydrostatic.hdf5 ]
+then
+    echo "Generating initial conditions for the isothermal potential box example..."
+    python makeIC.py 100000
+fi

 # Run for 10 dynamical times
-../swift -g -s -t 2 hydrostatic.yml 2>&1 | tee output.log
+../swift -g -s -t 1 hydrostatic.yml 2>&1 | tee output.log

 echo "Plotting density profiles"
 mkdir plots

--- a/examples/HydrostaticHalo/test_energy_conservation.py
+++ b/examples/HydrostaticHalo/test_energy_conservation.py
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
 import numpy as np
 import h5py as h5
-import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
 import sys

 n_snaps = int(sys.argv[1])
@@ -24,7 +45,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
@@ -45,7 +66,8 @@ for i in range(n_snaps):
    f = h5.File(filename,'r')
    coords_dset = f["PartType0/Coordinates"]
    coords = np.array(coords_dset)
-#translate coords by centre of box
+
+    #translate coords by centre of box
    header = f["Header"]
    snap_time = header.attrs["Time"]
    snap_time_cgs = snap_time * unit_time_cgs
@@ -73,7 +95,6 @@ for i in range(n_snaps):
    internal_energy_array = np.append(internal_energy_array,total_internal_energy)

 #put energies in units of v_c^2 and rescale by number of particles
-
 pe = potential_energy_array / (N*v_c**2)
 ke = kinetic_energy_array / (N*v_c**2)
 ie = internal_energy_array / (N*v_c**2)
@@ -82,14 +103,15 @@ te = pe + ke + ie
 dyn_time_cgs = r_vir_cgs / v_c_cgs
 time_array = time_array_cgs / dyn_time_cgs

-plt.plot(time_array,ke,label = "Kinetic Energy")
-plt.plot(time_array,pe,label = "Potential Energy")
-plt.plot(time_array,ie,label = "Internal Energy")
-plt.plot(time_array,te,label = "Total Energy")
-plt.legend(loc = "lower right")
-plt.xlabel(r"$t / t_{dyn}$")
-plt.ylabel(r"$E / v_c^2$")
-plt.title(r"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(N,v_c))
-plt.ylim((-2,2))
-plt.savefig("energy_conservation.png",format = 'png')
+figure()
+plot(time_array,ke,label = "Kinetic Energy")
+plot(time_array,pe,label = "Potential Energy")
+plot(time_array,ie,label = "Internal Energy")
+plot(time_array,te,label = "Total Energy")
+legend(loc = "lower right")
+xlabel(r"$t / t_{dyn}$")
+ylabel(r"$E / v_c^2$")
+title(r"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(N,v_c))
+ylim((-2,2))
+savefig("energy_conservation.png",format = 'png')

--- a/examples/HydrostaticHalo/velocity_profile.py
+++ b/examples/HydrostaticHalo/velocity_profile.py
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
 import numpy as np
 import h5py as h5
-import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
 import sys

 def do_binning(x,y,x_bin_edges):
@@ -62,7 +83,8 @@ for i in range(n_snaps):
    f = h5.File(filename,'r')
    coords_dset = f["PartType0/Coordinates"]
    coords = np.array(coords_dset)
-#translate coords by centre of box
+
+    #translate coords by centre of box
    header = f["Header"]
    snap_time = header.attrs["Time"]
    snap_time_cgs = snap_time * unit_time_cgs
@@ -73,16 +95,15 @@ for i in range(n_snaps):
    radius_cgs = radius*unit_length_cgs
    radius_over_virial_radius = radius_cgs / r_vir_cgs

-#get the internal energies
+    #get the internal energies
    vel_dset = f["PartType0/Velocities"]
    vel = np.array(vel_dset)

-#make dimensionless
+    #make dimensionless
    vel /= v_c
    r = radius_over_virial_radius

    #find radial component of velocity
-
    v_r = np.zeros(r.size)
    for j in range(r.size):
        v_r[j] = -np.dot(coords[j,:],vel[j,:])/radius[j]
@@ -94,18 +115,13 @@ for i in range(n_snaps):
    radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins) 
    binned_v_r = v_r_totals / hist

-    #calculate cooling radius
-
-    #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
-
-    plt.plot(radial_bin_mids,binned_v_r,'ko',label = "Average radial velocity in shell")
-    #plt.plot((0,1),(1,1),label = "Analytic Solution")
-    #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,2),'r',label = "Cooling radius")
-    plt.legend(loc = "upper right")
-    plt.xlabel(r"$r / r_{vir}$")
-    plt.ylabel(r"$v_r / v_c$")
-    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    plt.ylim((0,2))
+    figure()
+    plot(radial_bin_mids,binned_v_r,'ko',label = "Average radial velocity in shell")
+    legend(loc = "upper right")
+    xlabel(r"$r / r_{vir}$")
+    ylabel(r"$v_r / v_c$")
+    title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
+    ylim((-1,1))
    plot_filename = "./plots/radial_velocity_profile/velocity_profile_%03d.png" %i
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
+    savefig(plot_filename,format = "png")
+    close()
--- a/examples/UniformDMBox/makeIC.py
+++ b/examples/UniformDMBox/makeIC.py
@@ -26,7 +26,7 @@ from numpy import *
 # with a density of 1

 # Parameters
-periodic= 0           # 1 For periodic box
+periodic= 1           # 1 For periodic box
 boxSize = 1.
 rho = 1.
 L = int(sys.argv[1])  # Number of particles along one axis

--- a/examples/UniformDMBox/uniformBox.yml
+++ b/examples/UniformDMBox/uniformBox.yml
@@ -35,4 +35,4 @@ Statistics:

 # Parameters related to the initial conditions
 InitialConditions:
-  file_name:  ./uniformDMBox_100.hdf5     # The file to read
+  file_name:  ./uniformDMBox_50.hdf5     # The file to read
--- a/examples/plot_tasks.py
+++ b/examples/plot_tasks.py
 #!/usr/bin/env python
 """
 Usage:
-    plot_tasks.py input.dat output.png [time-range-ms]
+    plot_tasks.py [options] input.dat output.png

-where input.dat is a thread info file for a step.  Use the '-y interval'
-flag of the swift MPI commands to create these. The output plot will be
-called 'output.png'. Use the time-range-ms in millisecs to produce
-plots with the same time span.
+where input.dat is a thread info file for a step.  Use the '-y interval' flag
+of the swift command to create these. The output plot will be called
+'output.png'. The --limit option can be used to produce plots with the same
+time span and the --expand option to expand each thread line into '*expand'
+lines, so that adjacent tasks of the same type can be distinguished. Other
+options can be seen using the --help flag.

 This file is part of SWIFT.
 Copyright (c) 2015 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
                   Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
                   Matthieu Schaller (matthieu.schaller@durham.ac.uk)
-          (c) 2016 Peter W. Draper (p.w.draper@durham.ac.uk)
+          (c) 2017 Peter W. Draper (p.w.draper@durham.ac.uk)

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published
@@ -29,11 +31,42 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """

 import matplotlib
+matplotlib.use("Agg")
 import matplotlib.collections as collections
-matplotlib.use('Agg')
+import matplotlib.ticker as plticker
 import pylab as pl
-import numpy as np
 import sys
+import argparse
+
+#  Handle the command line.
+parser = argparse.ArgumentParser(description="Plot task graphs")
+
+parser.add_argument("input", help="Thread data file (-y output)")
+parser.add_argument("outpng", help="Name for output graphic file (PNG)")
+parser.add_argument("-l", "--limit", dest="limit",
+                    help="Upper time limit in millisecs (def: depends on data)",
+                    default=0, type=int)
+parser.add_argument("-e", "--expand", dest="expand",
+                    help="Thread expansion factor (def: 1)",
+                    default=1, type=int)
+parser.add_argument("--height", dest="height",
+                    help="Height of plot in inches (def: 4)",
+                    default=4., type=float)
+parser.add_argument("--width", dest="width",
+                    help="Width of plot in inches (def: 16)",
+                    default=16., type=float)
+parser.add_argument("--nolegend", dest="nolegend",
+                    help="Whether to show the legend (def: False)",
+                    default=False, action="store_true")
+parser.add_argument("-v", "--verbose", dest="verbose",
+                    help="Show colour assignments and other details (def: False)",
+                    default=False, action="store_true")
+
+args = parser.parse_args()
+infile = args.input
+outpng = args.outpng
+delta_t = args.limit
+expand = args.expand

 #  Basic plot configuration.
 PLOT_PARAMS = {"axes.labelsize": 10,
@@ -42,7 +75,7 @@ PLOT_PARAMS = {"axes.labelsize": 10,
               "legend.fontsize": 12,
               "xtick.labelsize": 10,
               "ytick.labelsize": 10,
-               "figure.figsize" : (16., 4.),
+               "figure.figsize" : (args.width, args.height),
               "figure.subplot.left" : 0.03,
               "figure.subplot.right" : 0.995,
               "figure.subplot.bottom" : 0.1,
@@ -56,9 +89,10 @@ pl.rcParams.update(PLOT_PARAMS)

 #  Tasks and subtypes. Indexed as in tasks.h.
 TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
-             "init", "ghost", "extra_ghost", "drift", "kick1", "kick2",
+             "init_grav", "ghost", "extra_ghost", "drift", "kick1", "kick2",
             "timestep", "send", "recv", "grav_top_level", "grav_long_range",
             "grav_mm", "grav_down", "cooling", "sourceterms", "count"]
+
 SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
            "tend", "xv", "rho", "gpart", "multipole", "spart", "count"]

@@ -69,14 +103,16 @@ FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
             "sub_pair/density", "recv/xv", "send/xv", "recv/rho", "send/rho",
             "recv/tend", "send/tend"]

-#  Get a number of colours for the various types.
-colours = ["black", "gray", "rosybrown", "firebrick", "red", "darksalmon",
-           "sienna", "sandybrown", "bisque", "tan", "moccasin", "gold", "darkkhaki",
-           "lightgoldenrodyellow", "olivedrab", "chartreuse", "darksage", "lightgreen",
-           "green", "mediumseagreen", "mediumaquamarine", "mediumturquoise", "darkslategrey",
-           "cyan", "cadetblue", "skyblue", "dodgerblue", "slategray", "darkblue",
-           "slateblue", "blueviolet", "mediumorchid", "purple", "magenta", "hotpink",
-           "pink"]
+#  A number of colours for the various types. Recycled when there are
+#  more task types than colours...
+
+colours = ["cyan", "lightgray", "darkblue", "yellow", "tan", "dodgerblue",
+           "sienna", "aquamarine", "bisque", "blue", "green", "brown",
+           "purple", "mocassin", "olivedrab", "chartreuse", "darksage",
+           "darkgreen", "green", "mediumseagreen", "mediumaquamarine",
+           "darkslategrey", "mediumturquoise", "black", "cadetblue", "skyblue",
+           "red", "slategray", "gold", "slateblue", "blueviolet",
+           "mediumorchid", "firebrick", "magenta", "hotpink", "pink"]
 maxcolours = len(colours)

 #  Set colours of task/subtype.
@@ -87,30 +123,21 @@ for task in TASKTYPES:
    ncolours = (ncolours + 1) % maxcolours

 SUBCOLOURS = {}
-for task in SUBTYPES:
+for task in FULLTYPES:
    SUBCOLOURS[task] = colours[ncolours]
    ncolours = (ncolours + 1) % maxcolours

-for task in FULLTYPES:
+for task in SUBTYPES:
    SUBCOLOURS[task] = colours[ncolours]
    ncolours = (ncolours + 1) % maxcolours

-#  Show docs if help is requested.
-if len( sys.argv ) == 2 and ( sys.argv[1][0:2] == "-h" or sys.argv[1][0:3] == "--h" ):
-    from pydoc import help
-    help( "__main__" )
-    sys.exit( 0 )
-
-#  Handle command-line.
-if len( sys.argv ) != 3 and len( sys.argv ) != 4:
-    print "Usage: ", sys.argv[0], "input.dat output.png [time-range-ms]"
-    sys.exit(1)
-
-infile = sys.argv[1]
-outpng = sys.argv[2]
-delta_t = 0
-if len( sys.argv ) == 4:
-    delta_t = int(sys.argv[3])
+#  For fiddling with colours...
+if args.verbose:
+    print "#Selected colours:"
+    for task in TASKCOLOURS.keys():
+        print "# " + task + ": " + TASKCOLOURS[task]
+    for task in SUBCOLOURS.keys():
+        print "# " + task + ": " + SUBCOLOURS[task]

 #  Read input.
 data = pl.loadtxt( infile )
@@ -118,51 +145,61 @@ data = pl.loadtxt( infile )
 nthread = int(max(data[:,0])) + 1
 print "Number of threads:", nthread

-# Recover the start and end time
+#  Recover the start and end time
 full_step = data[0,:]
 tic_step = int(full_step[4])