examples/GravityTests/DiscPatch/GravityOnly/test.pro

-;
-;  test energy / angular momentum conservation of test problem
-;
-
-iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare final and initial energy
-
-; set physical constants
-@physunits
-
-indir    = './'
-basefile = 'Disc-Patch_'
-
-; set properties of potential
-uL   = phys.pc                  ; unit of length
-uM   = phys.msun                ; unit of mass
-uV   = 1d5                      ; unit of velocity
-
-; properties of patch
-surface_density = 10.
-scale_height    = 100.
-
-; derived units
-constG   = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
-pcentre  = [0.,0.,300.] * pc / uL
-
-;
-infile = indir + basefile + '*'
-spawn,'ls -1 '+infile,res
-nfiles = n_elements(res)
-
-
-; choose: calculate change of energy and Lz, comparing first and last
-; snapshots for all particles, or do so for a subset
-
-; compare all
-ifile   = 0
-inf     = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
-id      = h5rd(inf,'PartType1/ParticleIDs')
-nfollow = n_elements(id)
-
-; follow a subset
-nfollow  = 500                    ; number of particles to follow
-
-;
-if (iplot eq 1) then begin
-   nskip = 1
-   nsave = nfiles
-endif else begin
-   nskip = nfiles - 2
-   nsave = 2
-endelse
-
-;
-lout     = fltarr(nfollow, nsave) ; Lz
-xout     = fltarr(nfollow, nsave) ; x
-yout     = fltarr(nfollow, nsave) ; y
-zout     = fltarr(nfollow, nsave) ; z
-eout     = fltarr(nfollow, nsave) ; energies
-ekin     = fltarr(nfollow, nsave)
-epot     = fltarr(nfollow, nsave) ; 2 pi G Sigma b ln(cosh(z/b)) + const
-tout     = fltarr(nsave)
-
-
-
-ifile  = 0
-isave = 0
-for ifile=0,nfiles-1,nskip do begin
-   inf    = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
-   time   = h5ra(inf, 'Header','Time')
-   p      = h5rd(inf,'PartType1/Coordinates')
-   v      = h5rd(inf,'PartType1/Velocities')
-   id     = h5rd(inf,'PartType1/ParticleIDs')
-   indx   = sort(id)
-
-;; ;  if you want to sort particles by ID
-;;    id     = id[indx]
-;;    for ic=0,2 do begin
-;;       tmp = reform(p[ic,*]) & p[ic,*] = tmp[indx]
-;;       tmp = reform(v[ic,*]) & v[ic,*] = tmp[indx]
-;;    endfor
-   
-
-; calculate energy
-   dd  = size(p,/dimen) & npart = dd[1]
-   ener = fltarr(npart)
-   dr   = fltarr(npart) & dv = dr
-   for ic=0,2 do dr[*] = dr[*] + (p[ic,*]-pcentre[ic])^2
-   for ic=0,2 do dv[*] = dv[*] + v[ic,*]^2
-   xout[*,isave] = p[0,0:nfollow-1]-pcentre[0]
-   yout[*,isave] = p[1,0:nfollow-1]-pcentre[1]
-   zout[*,isave] = p[2,0:nfollow-1]-pcentre[2]
-   Lz  = (p[0,*]-pcentre[0]) * v[1,*] - (p[1,*]-pcentre[1]) * v[0,*]
-   dz  = reform(p[2,0:nfollow-1]-pcentre[2])
-;   print,'time = ',time,p[0,0],v[0,0],id[0]
-   ek   = 0.5 * dv
-   ep   = fltarr(nfollow)
-   ep   = 2 * !pi * constG * surface_density * scale_height * alog(cosh(abs(dz)/scale_height))
-   ener = ek + ep
-   tout(isave) = time
-   lout[*,isave] = lz[0:nfollow-1]
-   eout(*,isave) = ener[0:nfollow-1]
-   ekin(*,isave) = ek[0:nfollow-1]
-   epot(*,isave) = ep[0:nfollow-1]
-   print,format='('' time= '',f7.1,'' E= '',f9.2,'' Lz= '',e9.2)', time,eout[0],lz[0]
-   isave = isave + 1
-   
-endfor
-
-x0 = reform(xout[0,*])
-y0 = reform(xout[1,*])
-z0 = reform(xout[2,*])
-
-; calculate relative energy change
-de    = 0.0 * eout
-dl    = 0.0 * lout
-nsave = isave
-for ifile=1, nsave-1 do de[*,ifile] = (eout[*,ifile]-eout[*,0])/eout[*,0]
-for ifile=1, nsave-1 do dl[*,ifile] = (lout[*,ifile] - lout[*,0])/lout[*,0]
-
-
-; calculate statistics of energy changes
-print,' relatve energy change: (per cent) ',minmax(de) * 100.
-print,' relative Lz    change: (per cent) ',minmax(dl) * 100.
-
-; plot enery and Lz conservation for some particles
-if(iplot eq 1) then begin
-; plot results on energy conservation for some particles
-   nplot = min(10, nfollow)
-   win,0
-   xr = [min(tout), max(tout)]
-   yr = [-2,2]*1d-2             ; in percent
-   plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='time',ytitle='dE/E, dL/L (%)'
-   for i=0,nplot-1 do oplot,tout,de[i,*]
-   for i=0,nplot-1 do oplot,tout,dl[i,*],color=red
-   legend,['dE/E','dL/L'],linestyle=[0,0],color=[black,red],box=0,/bottom,/left
-   screen_to_png,'e-time.png'
-
-;  plot vertical oscillation
-   win,2
-   xr = [min(tout), max(tout)]
-   yr = [-3,3]*scale_height
-   plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/iso,/nodata,xtitle='x',ytitle='y'
-   color = floor(findgen(nplot)*255/float(nplot))
-   for i=0,nplot-1 do oplot,tout,zout[i,*],color=color(i)
-   screen_to_png,'orbit.png'
-
-; make histogram of energy changes at end
-   win,6
-   ohist,de,x,y,-0.05,0.05,0.001
-   plot,x,y,psym=10,xtitle='de (%)'
-   screen_to_png,'de-hist.png'
-
-
-endif
-
-end
-
-

examples/GravityTests/DiscPatch/GravityOnly/test.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# 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 h5py as h5
+import numpy as np
+import math
+
+num_snapshots = 61
+
+# physical constants in cgs
+NEWTON_GRAVITY_CGS = 6.67430e-8
+SOLAR_MASS_IN_CGS = 1.98841e33
+PARSEC_IN_CGS = 3.08567758149e18
+YEAR_IN_CGS = 3.15569251e7
+
+# choice of units
+const_unit_length_in_cgs = PARSEC_IN_CGS
+const_unit_mass_in_cgs = SOLAR_MASS_IN_CGS
+const_unit_velocity_in_cgs = 1e5
+
+# parameters of potential
+surface_density = 10.0
+scale_height = 100.0
+centre = np.array([300.0, 300.0, 300.0])
+
+# constants
+const_unit_time_in_cgs = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+const_G = (
+    NEWTON_GRAVITY_CGS
+    * const_unit_mass_in_cgs
+    * const_unit_time_in_cgs
+    * const_unit_time_in_cgs
+    / (const_unit_length_in_cgs * const_unit_length_in_cgs * const_unit_length_in_cgs)
+)
+
+
+t = np.zeros(num_snapshots)
+E_k = np.zeros(num_snapshots)
+E_p = np.zeros(num_snapshots)
+E_tot = np.zeros(num_snapshots)
+
+# loop over the snapshots
+for i in range(num_snapshots):
+
+    filename = "Disc-Patch_%04d.hdf5" % i
+    f = h5.File(filename, "r")
+
+    # time
+    t[i] = f["/Header"].attrs.get("Time")[0]
+
+    # positions and velocities of the particles
+    p = f["/PartType1/Coordinates"][:]
+    v = f["/PartType1/Velocities"][:]
+
+    # Kinetic energy
+    E_k[i] = np.sum(0.5 * (v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2))
+
+    # Potential energy
+    d = p[:, 0] - centre[0]
+    E_p[i] = np.sum(
+        2.0
+        * math.pi
+        * const_G
+        * surface_density
+        * scale_height
+        * np.log(np.cosh(np.abs(d) / scale_height))
+    )
+
+    # Total energy
+    E_tot[i] = E_k[i] + E_p[i]
+
+
+# Maximal change
+delta_E = np.max(np.abs(E_tot - E_tot[0])) / E_tot[0]
+
+print(("Maximal relative energy change   :", delta_E * 100, "%"))

examples/GravityTests/DiscPatch/HydroStatic/getGlass.sh

 #!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5

examples/GravityTests/DiscPatch/HydroStatic/makeIC.py

@@ -2,7 +2,7 @@
 # This file is part of SWIFT.
 # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
 #                    Tom Theuns (tom.theuns@durham.ac.uk)
-#               2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#               2017 Matthieu Schaller (schaller@strw.leidenuniv.nl)
 #                    Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
 #
 # This program is free software: you can redistribute it and/or modify
@@ -39,13 +39,13 @@ import random
 #                         Size of the patch -- side_length
 
 # Parameters of the gas disc
-surface_density = 10.
-scale_height    = 100.
-gas_gamma       = 5./3.
+surface_density = 10.0
+scale_height = 100.0
+gas_gamma = 5.0 / 3.0
 
 # Parameters of the problem
-x_factor        = 2
-side_length     = 400.
+x_factor = 2
+side_length = 400.0
 
 # File
 fileName = "Disc-Patch.hdf5"
@@ -53,55 +53,62 @@ fileName = "Disc-Patch.hdf5"
 ####################################################################
 
 # physical constants in cgs
-NEWTON_GRAVITY_CGS  = 6.67408e-8
-SOLAR_MASS_IN_CGS   = 1.98848e33
-PARSEC_IN_CGS       = 3.08567758e18
-PROTON_MASS_IN_CGS  = 1.672621898e-24
-BOLTZMANN_IN_CGS    = 1.38064852e-16
-YEAR_IN_CGS         = 3.15569252e7
+NEWTON_GRAVITY_CGS = 6.67408e-8
+SOLAR_MASS_IN_CGS = 1.98848e33
+PARSEC_IN_CGS = 3.08567758e18
+PROTON_MASS_IN_CGS = 1.672621898e-24
+BOLTZMANN_IN_CGS = 1.38064852e-16
+YEAR_IN_CGS = 3.15569252e7
 
 # choice of units
-unit_length_in_cgs   =   (PARSEC_IN_CGS)
-unit_mass_in_cgs     =   (SOLAR_MASS_IN_CGS)
-unit_velocity_in_cgs =   (1e5)
-unit_time_in_cgs     =   unit_length_in_cgs / unit_velocity_in_cgs
+unit_length_in_cgs = PARSEC_IN_CGS
+unit_mass_in_cgs = SOLAR_MASS_IN_CGS
+unit_velocity_in_cgs = 1e5
+unit_time_in_cgs = unit_length_in_cgs / unit_velocity_in_cgs
 
-print "UnitMass_in_cgs:     %.5e"%unit_mass_in_cgs
-print "UnitLength_in_cgs:   %.5e"%unit_length_in_cgs
-print "UnitVelocity_in_cgs: %.5e"%unit_velocity_in_cgs
-print "UnitTime_in_cgs:     %.5e"%unit_time_in_cgs
-print ""
+print("UnitMass_in_cgs:     %.5e" % unit_mass_in_cgs)
+print("UnitLength_in_cgs:   %.5e" % unit_length_in_cgs)
+print("UnitVelocity_in_cgs: %.5e" % unit_velocity_in_cgs)
+print("UnitTime_in_cgs:     %.5e" % unit_time_in_cgs)
+print("")
 
 # Derived units
-const_G  = NEWTON_GRAVITY_CGS * unit_mass_in_cgs * unit_time_in_cgs**2 * \
-           unit_length_in_cgs**-3
-const_mp = PROTON_MASS_IN_CGS * unit_mass_in_cgs**-1
-const_kb = BOLTZMANN_IN_CGS * unit_mass_in_cgs**-1 * unit_length_in_cgs**-2 * \
-           unit_time_in_cgs**2
-
-print "--- Some constants [internal units] ---"
-print "G_Newton:    %.5e"%const_G
-print "m_proton:    %.5e"%const_mp
-print "k_boltzmann: %.5e"%const_kb
-print ""
+const_G = (
+    NEWTON_GRAVITY_CGS
+    * unit_mass_in_cgs
+    * unit_time_in_cgs ** 2
+    * unit_length_in_cgs ** -3
+)
+const_mp = PROTON_MASS_IN_CGS * unit_mass_in_cgs ** -1
+const_kb = (
+    BOLTZMANN_IN_CGS
+    * unit_mass_in_cgs ** -1
+    * unit_length_in_cgs ** -2
+    * unit_time_in_cgs ** 2
+)
+
+print("--- Some constants [internal units] ---")
+print("G_Newton:    %.5e" % const_G)
+print("m_proton:    %.5e" % const_mp)
+print("k_boltzmann: %.5e" % const_kb)
+print("")
 
 # derived quantities
-temp       = math.pi * const_G * surface_density * scale_height * const_mp / \
-             const_kb
-u_therm    = const_kb * temp / ((gas_gamma-1) * const_mp)
-v_disp     = math.sqrt(2 * u_therm)
-soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma-1.)))
-t_dyn      = math.sqrt(scale_height / (const_G * surface_density))
-t_cross    = scale_height / soundspeed
-
-print "--- Properties of the gas [internal units] ---"
-print "Gas temperature:     %.5e"%temp
-print "Gas thermal_energy:  %.5e"%u_therm
-print "Dynamical time:      %.5e"%t_dyn
-print "Sound crossing time: %.5e"%t_cross
-print "Gas sound speed:     %.5e"%soundspeed
-print "Gas 3D vel_disp:     %.5e"%v_disp
-print ""
+temp = math.pi * const_G * surface_density * scale_height * const_mp / const_kb
+u_therm = const_kb * temp / ((gas_gamma - 1) * const_mp)
+v_disp = math.sqrt(2 * u_therm)
+soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma - 1.0)))
+t_dyn = math.sqrt(scale_height / (const_G * surface_density))
+t_cross = scale_height / soundspeed
+
+print("--- Properties of the gas [internal units] ---")
+print("Gas temperature:     %.5e" % temp)
+print("Gas thermal_energy:  %.5e" % u_therm)
+print("Dynamical time:      %.5e" % t_dyn)
+print("Sound crossing time: %.5e" % t_cross)
+print("Gas sound speed:     %.5e" % soundspeed)
+print("Gas 3D vel_disp:     %.5e" % v_disp)
+print("")
 
 # Problem properties
 boxSize_x = side_length
@@ -109,71 +116,75 @@ boxSize_y = boxSize_x
 boxSize_z = boxSize_x
 boxSize_x *= x_factor
 volume = boxSize_x * boxSize_y * boxSize_z
-M_tot = boxSize_y * boxSize_z * surface_density * \
-        math.tanh(boxSize_x / (2. * scale_height))
+M_tot = (
+    boxSize_y
+    * boxSize_z
+    * surface_density
+    * math.tanh(boxSize_x / (2.0 * scale_height))
+)
 density = M_tot / volume
-entropy = (gas_gamma - 1.) * u_therm / density**(gas_gamma - 1.)
+entropy = (gas_gamma - 1.0) * u_therm / density ** (gas_gamma - 1.0)
 
-print "--- Problem properties [internal units] ---"
-print "Box:        [%.1f, %.1f, %.1f]"%(boxSize_x, boxSize_y, boxSize_z)
-print "Volume:     %.5e"%volume
-print "Total mass: %.5e"%M_tot
-print "Density:    %.5e"%density
-print "Entropy:    %.5e"%entropy
-print ""
+print("--- Problem properties [internal units] ---")
+print("Box:        [%.1f, %.1f, %.1f]" % (boxSize_x, boxSize_y, boxSize_z))
+print("Volume:     %.5e" % volume)
+print("Total mass: %.5e" % M_tot)
+print("Density:    %.5e" % density)
+print("Entropy:    %.5e" % entropy)
+print("")
 
 ####################################################################
 
 # Read glass pre-ICs
-infile  = h5py.File('glassCube_32.hdf5', "r")
-one_glass_pos = infile["/PartType0/Coordinates"][:,:]
-one_glass_h   = infile["/PartType0/SmoothingLength"][:]
+infile = h5py.File("glassCube_32.hdf5", "r")
+one_glass_pos = infile["/PartType0/Coordinates"][:, :]
+one_glass_h = infile["/PartType0/SmoothingLength"][:]
 
 # Rescale to the problem size
 one_glass_pos *= side_length
-one_glass_h   *= side_length
+one_glass_h *= side_length
 
 # Now create enough copies to fill the volume in x
 pos = np.copy(one_glass_pos)
 h = np.copy(one_glass_h)
 for i in range(1, x_factor):
-    one_glass_pos[:,0] += side_length
+    one_glass_pos[:, 0] += side_length
     pos = np.append(pos, one_glass_pos, axis=0)
-    h   = np.append(h, one_glass_h, axis=0)
+    h = np.append(h, one_glass_h, axis=0)
 
 # Compute further properties of ICs
 numPart = np.size(h)
 mass = M_tot / numPart
 
-print "--- Particle properties [internal units] ---"
-print "Number part.: ", numPart
-print "Part. mass:   %.5e"%mass
-print ""
+print("--- Particle properties [internal units] ---")
+print("Number part.: ", numPart)
+print("Part. mass:   %.5e" % mass)
+print("")
 
 # Create additional arrays
-u    = np.ones(numPart) * u_therm
+u = np.ones(numPart) * u_therm
 mass = np.ones(numPart) * mass
-vel  = np.zeros((numPart, 3))
-ids  = 1 + np.linspace(0, numPart, numPart, endpoint=False)
+vel = np.zeros((numPart, 3))
+ids = 1 + np.linspace(0, numPart, numPart, endpoint=False)
 
 ####################################################################
 # Create and write output file
 
-#File
-file = h5py.File(fileName, 'w')
+# File
+file = h5py.File(fileName, "w")
 
-#Units
+# Units
 grp = file.create_group("/Units")
 grp.attrs["Unit length in cgs (U_L)"] = unit_length_in_cgs
 grp.attrs["Unit mass in cgs (U_M)"] = unit_mass_in_cgs
 grp.attrs["Unit time in cgs (U_t)"] = unit_time_in_cgs
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
 # Header
 grp = file.create_group("/Header")
 grp.attrs["BoxSize"] = [boxSize_x, boxSize_y, boxSize_z]
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["Time"] = 0.0
@@ -183,20 +194,20 @@ grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["Dimension"] = 3
 
 # write gas particles
-grp0   = file.create_group("/PartType0")
+grp0 = file.create_group("/PartType0")
 
-ds = grp0.create_dataset('Coordinates', (numPart, 3), 'f', data=pos)
-ds = grp0.create_dataset('Velocities', (numPart, 3), 'f')
-ds = grp0.create_dataset('Masses', (numPart,), 'f', data=mass)
-ds = grp0.create_dataset('SmoothingLength', (numPart,), 'f', data=h)
-ds = grp0.create_dataset('InternalEnergy', (numPart,), 'f', data=u)
-ds = grp0.create_dataset('ParticleIDs', (numPart, ), 'L', data=ids)
+ds = grp0.create_dataset("Coordinates", (numPart, 3), "f", data=pos)
+ds = grp0.create_dataset("Velocities", (numPart, 3), "f")
+ds = grp0.create_dataset("Masses", (numPart,), "f", data=mass)
+ds = grp0.create_dataset("SmoothingLength", (numPart,), "f", data=h)
+ds = grp0.create_dataset("InternalEnergy", (numPart,), "f", data=u)
+ds = grp0.create_dataset("ParticleIDs", (numPart,), "L", data=ids)
 
 ####################################################################
 
-print "--- Runtime parameters (YAML file): ---"
-print "DiscPatchPotential:surface_density:    ", surface_density
-print "DiscPatchPotential:scale_height:       ", scale_height
-print "DiscPatchPotential:x_disc:             ", 0.5 * boxSize_x
-print "EoS:isothermal_internal_energy: %ef"%u_therm
-print ""
+print("--- Runtime parameters (YAML file): ---")
+print("DiscPatchPotential:surface_density:    ", surface_density)
+print("DiscPatchPotential:scale_height:       ", scale_height)
+print("DiscPatchPotential:x_disc:             ", 0.5 * boxSize_x)
+print("EoS:isothermal_internal_energy: %ef" % u_therm)
+print("")

examples/GravityTests/DiscPatch/HydroStatic/plotSolution.py

 ################################################################################
 # This file is part of SWIFT.
 # Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
-#                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#                    Matthieu Schaller (schaller@strw.leidenuniv.nl)
 #
 # 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
@@ -27,95 +27,98 @@
 import numpy as np
 import h5py
 import matplotlib
+
 matplotlib.use("Agg")
 import pylab as pl
 import glob
 import sys
 
 # Parameters
-surface_density = 10.
-scale_height = 100.
-x_disc = 400.
-x_trunc = 300.
-x_max = 350.
+surface_density = 10.0
+scale_height = 100.0
+x_disc = 400.0
+x_trunc = 300.0
+x_max = 350.0
 utherm = 20.2678457288
-gamma = 5. / 3.
+gamma = 5.0 / 3.0
 
 start = 0
 stop = 21
 if len(sys.argv) > 1:
-  start = int(sys.argv[1])
+    start = int(sys.argv[1])
 if len(sys.argv) > 2:
-  stop = int(sys.argv[2])
+    stop = int(sys.argv[2])
 
 # Get the analytic solution for the density
 def get_analytic_density(x):
-  return 0.5 * surface_density / scale_height / \
-           np.cosh( (x - x_disc) / scale_height )**2
+    return (
+        0.5 * surface_density / scale_height / np.cosh((x - x_disc) / scale_height) ** 2
+    )
+
 
 # Get the analytic solution for the (isothermal) pressure
 def get_analytic_pressure(x):
-  return (gamma - 1.) * utherm * get_analytic_density(x)
+    return (gamma - 1.0) * utherm * get_analytic_density(x)
+
 
 # Get the data fields to plot from the snapshot file with the given name:
 #  snapshot time, x-coord, density, pressure, velocity norm
 def get_data(name):
-  file = h5py.File(name, "r")
-  coords = np.array(file["/PartType0/Coordinates"])
-  rho = np.array(file["/PartType0/Density"])
-  u = np.array(file["/PartType0/InternalEnergy"])
-  v = np.array(file["/PartType0/Velocities"])
+    file = h5py.File(name, "r")
+    coords = np.array(file["/PartType0/Coordinates"])
+    rho = np.array(file["/PartType0/Densities"])
+    v = np.array(file["/PartType0/Velocities"])
+    P = np.array(file["/PartType0/Pressures"])
 
-  P = (gamma - 1.) * rho * u
+    vtot = np.sqrt(v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2)
 
-  vtot = np.sqrt( v[:,0]**2 + v[:,1]**2 + v[:,2]**2 )
+    return float(file["/Header"].attrs["Time"]), coords[:, 0], rho, P, vtot
 
-  return float(file["/Header"].attrs["Time"]), coords[:,0], rho, P, vtot
 
 # scan the folder for snapshot files and plot all of them (within the requested
 # range)
 for f in sorted(glob.glob("Disc-Patch_*.hdf5")):
-  num = int(f[-8:-5])
-  if num < start or num > stop:
-    continue
-
-  print "processing", f, "..."
-
-  xrange = np.linspace(0., 2. * x_disc, 1000)
-  time, x, rho, P, v = get_data(f)
-
-  fig, ax = pl.subplots(3, 1, sharex = True)
-
-  ax[0].plot(x, rho, "r.")
-  ax[0].plot(xrange, get_analytic_density(xrange), "k-")
-  ax[0].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
-  ax[0].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
-  ax[0].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[0].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[0].set_ylim(0., 1.2 * get_analytic_density(x_disc))
-  ax[0].set_ylabel("density")
-
-  ax[1].plot(x, v, "r.")
-  ax[1].plot(xrange, np.zeros(len(xrange)), "k-")
-  ax[1].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
-  ax[1].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
-  ax[1].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[1].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[1].set_ylim(-0.5, 10.)
-  ax[1].set_ylabel("velocity norm")
-
-  ax[2].plot(x, P, "r.")
-  ax[2].plot(xrange, get_analytic_pressure(xrange), "k-")
-  ax[2].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
-  ax[2].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
-  ax[2].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[2].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[2].set_xlim(0., 2. * x_disc)
-  ax[2].set_ylim(0., 1.2 * get_analytic_pressure(x_disc))
-  ax[2].set_xlabel("x")
-  ax[2].set_ylabel("pressure")
-
-  pl.suptitle("t = {0:.2f}".format(time))
-
-  pl.savefig("{name}.png".format(name = f[:-5]))
-  pl.close()
+    num = int(f[-8:-5])
+    if num < start or num > stop:
+        continue
+
+    print(("processing", f, "..."))
+
+    xrange = np.linspace(0.0, 2.0 * x_disc, 1000)
+    time, x, rho, P, v = get_data(f)
+
+    fig, ax = pl.subplots(3, 1, sharex=True)
+
+    ax[0].plot(x, rho, "r.")
+    ax[0].plot(xrange, get_analytic_density(xrange), "k-")
+    ax[0].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+    ax[0].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+    ax[0].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[0].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[0].set_ylim(0.0, 1.2 * get_analytic_density(x_disc))
+    ax[0].set_ylabel("density")
+
+    ax[1].plot(x, v, "r.")
+    ax[1].plot(xrange, np.zeros(len(xrange)), "k-")
+    ax[1].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+    ax[1].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+    ax[1].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[1].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[1].set_ylim(-0.5, 10.0)
+    ax[1].set_ylabel("velocity norm")
+
+    ax[2].plot(x, P, "r.")
+    ax[2].plot(xrange, get_analytic_pressure(xrange), "k-")
+    ax[2].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+    ax[2].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+    ax[2].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[2].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[2].set_xlim(0.0, 2.0 * x_disc)
+    ax[2].set_ylim(0.0, 1.2 * get_analytic_pressure(x_disc))
+    ax[2].set_xlabel("x")
+    ax[2].set_ylabel("pressure")
+
+    pl.suptitle("t = {0:.2f}".format(time))
+
+    pl.savefig("{name}.png".format(name=f[:-5]))
+    pl.close()

examples/GravityTests/DiscPatch/HydroStatic/run.sh

@@ -9,10 +9,10 @@ fi
 if [ ! -e Disc-Patch.hdf5 ]
 then
     echo "Generating initial conditions for the disc patch example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../../swift --external-gravity --hydro --threads=4 disc-patch-icc.yml 2>&1 | tee output.log
+../../../../swift --external-gravity --hydro --threads=4 disc-patch-icc.yml 2>&1 | tee output.log
 
-python plotSolution.py
+python3 plotSolution.py

examples/GravityTests/DiscPatch/HydroStatic_1D/makeIC.py

@@ -2,7 +2,7 @@
 # This file is part of SWIFT.
 # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
 #                    Tom Theuns (tom.theuns@durham.ac.uk)
-#               2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#               2017 Matthieu Schaller (schaller@strw.leidenuniv.nl)
 #                    Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
 #
 # This program is free software: you can redistribute it and/or modify
@@ -39,14 +39,14 @@ import random
 #                         Size of the patch -- side_length
 
 # Parameters of the gas disc
-surface_density = 10.
-scale_height    = 100.
-gas_gamma       = 5./3.
+surface_density = 10.0
+scale_height = 100.0
+gas_gamma = 5.0 / 3.0
 
 # Parameters of the problem
-x_factor        = 2
-side_length     = 400.
-numPart         = 1000
+x_factor = 2
+side_length = 400.0
+numPart = 1000
 
 # File
 fileName = "Disc-Patch.hdf5"
@@ -54,112 +54,119 @@ fileName = "Disc-Patch.hdf5"
 ####################################################################
 
 # physical constants in cgs
-NEWTON_GRAVITY_CGS  = 6.67408e-8
-SOLAR_MASS_IN_CGS   = 1.98848e33
-PARSEC_IN_CGS       = 3.08567758e18
-PROTON_MASS_IN_CGS  = 1.672621898e-24
-BOLTZMANN_IN_CGS    = 1.38064852e-16
-YEAR_IN_CGS         = 3.15569252e7
+NEWTON_GRAVITY_CGS = 6.67408e-8
+SOLAR_MASS_IN_CGS = 1.98848e33
+PARSEC_IN_CGS = 3.08567758e18
+PROTON_MASS_IN_CGS = 1.672621898e-24
+BOLTZMANN_IN_CGS = 1.38064852e-16
+YEAR_IN_CGS = 3.15569252e7
 
 # choice of units
-unit_length_in_cgs   =   (PARSEC_IN_CGS)
-unit_mass_in_cgs     =   (SOLAR_MASS_IN_CGS)
-unit_velocity_in_cgs =   (1e5)
-unit_time_in_cgs     =   unit_length_in_cgs / unit_velocity_in_cgs
+unit_length_in_cgs = PARSEC_IN_CGS
+unit_mass_in_cgs = SOLAR_MASS_IN_CGS
+unit_velocity_in_cgs = 1e5
+unit_time_in_cgs = unit_length_in_cgs / unit_velocity_in_cgs
 
-print "UnitMass_in_cgs:     %.5e"%unit_mass_in_cgs
-print "UnitLength_in_cgs:   %.5e"%unit_length_in_cgs
-print "UnitVelocity_in_cgs: %.5e"%unit_velocity_in_cgs
-print "UnitTime_in_cgs:     %.5e"%unit_time_in_cgs
-print ""
+print("UnitMass_in_cgs:     %.5e" % unit_mass_in_cgs)
+print("UnitLength_in_cgs:   %.5e" % unit_length_in_cgs)
+print("UnitVelocity_in_cgs: %.5e" % unit_velocity_in_cgs)
+print("UnitTime_in_cgs:     %.5e" % unit_time_in_cgs)
+print("")
 
 # Derived units
-const_G  = NEWTON_GRAVITY_CGS * unit_mass_in_cgs * unit_time_in_cgs**2 * \
-           unit_length_in_cgs**-3
-const_mp = PROTON_MASS_IN_CGS * unit_mass_in_cgs**-1
-const_kb = BOLTZMANN_IN_CGS * unit_mass_in_cgs**-1 * unit_length_in_cgs**-2 * \
-           unit_time_in_cgs**2
-
-print "--- Some constants [internal units] ---"
-print "G_Newton:    %.5e"%const_G
-print "m_proton:    %.5e"%const_mp
-print "k_boltzmann: %.5e"%const_kb
-print ""
+const_G = (
+    NEWTON_GRAVITY_CGS
+    * unit_mass_in_cgs
+    * unit_time_in_cgs ** 2
+    * unit_length_in_cgs ** -3
+)
+const_mp = PROTON_MASS_IN_CGS * unit_mass_in_cgs ** -1
+const_kb = (
+    BOLTZMANN_IN_CGS
+    * unit_mass_in_cgs ** -1
+    * unit_length_in_cgs ** -2
+    * unit_time_in_cgs ** 2
+)
+
+print("--- Some constants [internal units] ---")
+print("G_Newton:    %.5e" % const_G)
+print("m_proton:    %.5e" % const_mp)
+print("k_boltzmann: %.5e" % const_kb)
+print("")
 
 # derived quantities
-temp       = math.pi * const_G * surface_density * scale_height * const_mp / \
-             const_kb
-u_therm    = const_kb * temp / ((gas_gamma-1) * const_mp)
-v_disp     = math.sqrt(2 * u_therm)
-soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma-1.)))
-t_dyn      = math.sqrt(scale_height / (const_G * surface_density))
-t_cross    = scale_height / soundspeed
-
-print "--- Properties of the gas [internal units] ---"
-print "Gas temperature:     %.5e"%temp
-print "Gas thermal_energy:  %.5e"%u_therm
-print "Dynamical time:      %.5e"%t_dyn
-print "Sound crossing time: %.5e"%t_cross
-print "Gas sound speed:     %.5e"%soundspeed
-print "Gas 3D vel_disp:     %.5e"%v_disp
-print ""
+temp = math.pi * const_G * surface_density * scale_height * const_mp / const_kb
+u_therm = const_kb * temp / ((gas_gamma - 1) * const_mp)
+v_disp = math.sqrt(2 * u_therm)
+soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma - 1.0)))
+t_dyn = math.sqrt(scale_height / (const_G * surface_density))
+t_cross = scale_height / soundspeed
+
+print("--- Properties of the gas [internal units] ---")
+print("Gas temperature:     %.5e" % temp)
+print("Gas thermal_energy:  %.5e" % u_therm)
+print("Dynamical time:      %.5e" % t_dyn)
+print("Sound crossing time: %.5e" % t_cross)
+print("Gas sound speed:     %.5e" % soundspeed)
+print("Gas 3D vel_disp:     %.5e" % v_disp)
+print("")
 
 # Problem properties
 boxSize_x = side_length
 boxSize_x *= x_factor
 volume = boxSize_x
-M_tot = surface_density * math.tanh(boxSize_x / (2. * scale_height))
+M_tot = surface_density * math.tanh(boxSize_x / (2.0 * scale_height))
 density = M_tot / volume
-entropy = (gas_gamma - 1.) * u_therm / density**(gas_gamma - 1.)
+entropy = (gas_gamma - 1.0) * u_therm / density ** (gas_gamma - 1.0)
 
-print "--- Problem properties [internal units] ---"
-print "Box:        %.1f"%boxSize_x
-print "Volume:     %.5e"%volume
-print "Total mass: %.5e"%M_tot
-print "Density:    %.5e"%density
-print "Entropy:    %.5e"%entropy
-print ""
+print("--- Problem properties [internal units] ---")
+print("Box:        %.1f" % boxSize_x)
+print("Volume:     %.5e" % volume)
+print("Total mass: %.5e" % M_tot)
+print("Density:    %.5e" % density)
+print("Entropy:    %.5e" % entropy)
+print("")
 
 ####################################################################
 
 # Now create enough copies to fill the volume in x
 pos = np.zeros((numPart, 3))
-h = np.zeros(numPart) + 2. * boxSize_x / numPart
+h = np.zeros(numPart) + 2.0 * boxSize_x / numPart
 for i in range(numPart):
     pos[i, 0] = (i + 0.5) * boxSize_x / numPart
 
 # Compute further properties of ICs
 mass = M_tot / numPart
 
-print "--- Particle properties [internal units] ---"
-print "Number part.: ", numPart
-print "Part. mass:   %.5e"%mass
-print ""
+print("--- Particle properties [internal units] ---")
+print("Number part.: ", numPart)
+print("Part. mass:   %.5e" % mass)
+print("")
 
 # Create additional arrays
-u    = np.ones(numPart) * u_therm
+u = np.ones(numPart) * u_therm
 mass = np.ones(numPart) * mass
-vel  = np.zeros((numPart, 3))
-ids  = 1 + np.linspace(0, numPart, numPart, endpoint=False)
+vel = np.zeros((numPart, 3))
+ids = 1 + np.linspace(0, numPart, numPart, endpoint=False)
 
 ####################################################################
 # Create and write output file
 
-#File
-file = h5py.File(fileName, 'w')
+# File
+file = h5py.File(fileName, "w")
 
-#Units
+# Units
 grp = file.create_group("/Units")
 grp.attrs["Unit length in cgs (U_L)"] = unit_length_in_cgs
 grp.attrs["Unit mass in cgs (U_M)"] = unit_mass_in_cgs
 grp.attrs["Unit time in cgs (U_t)"] = unit_time_in_cgs
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
 # Header
 grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = [boxSize_x, 1., 1.]
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["BoxSize"] = [boxSize_x, 1.0, 1.0]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["Time"] = 0.0
@@ -169,22 +176,22 @@ grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["Dimension"] = 1
 
 # write gas particles
-grp0   = file.create_group("/PartType0")
+grp0 = file.create_group("/PartType0")
 
-ds = grp0.create_dataset('Coordinates', (numPart, 3), 'f', data=pos)
-ds = grp0.create_dataset('Velocities', (numPart, 3), 'f')
-ds = grp0.create_dataset('Masses', (numPart,), 'f', data=mass)
-ds = grp0.create_dataset('SmoothingLength', (numPart,), 'f', data=h)
-ds = grp0.create_dataset('InternalEnergy', (numPart,), 'f', data=u)
-ds = grp0.create_dataset('ParticleIDs', (numPart, ), 'L', data=ids)
+ds = grp0.create_dataset("Coordinates", (numPart, 3), "f", data=pos)
+ds = grp0.create_dataset("Velocities", (numPart, 3), "f")
+ds = grp0.create_dataset("Masses", (numPart,), "f", data=mass)
+ds = grp0.create_dataset("SmoothingLength", (numPart,), "f", data=h)
+ds = grp0.create_dataset("InternalEnergy", (numPart,), "f", data=u)
+ds = grp0.create_dataset("ParticleIDs", (numPart,), "L", data=ids)
 
 ####################################################################
 
-print "--- Runtime parameters (YAML file): ---"
-print "DiscPatchPotential:surface_density:    ", surface_density
-print "DiscPatchPotential:scale_height:       ", scale_height
-print "DiscPatchPotential:x_disc:             ", 0.5 * boxSize_x
-print ""
+print("--- Runtime parameters (YAML file): ---")
+print("DiscPatchPotential:surface_density:    ", surface_density)
+print("DiscPatchPotential:scale_height:       ", scale_height)
+print("DiscPatchPotential:x_disc:             ", 0.5 * boxSize_x)
+print("")
 
-print "--- Constant parameters: ---"
-print "const_isothermal_internal_energy: %ef"%u_therm
+print("--- Constant parameters: ---")
+print("const_isothermal_internal_energy: %ef" % u_therm)

examples/GravityTests/DiscPatch/HydroStatic_1D/plotSolution.py

 ################################################################################
 # This file is part of SWIFT.
 # Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
-#                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#                    Matthieu Schaller (schaller@strw.leidenuniv.nl)
 #
 # 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
@@ -27,95 +27,100 @@
 import numpy as np
 import h5py
 import matplotlib
+
 matplotlib.use("Agg")
 import pylab as pl
 import glob
 import sys
 
 # Parameters
-surface_density = 10.
-scale_height = 100.
-x_disc = 400.
-x_trunc = 300.
-x_max = 350.
+surface_density = 10.0
+scale_height = 100.0
+x_disc = 400.0
+x_trunc = 300.0
+x_max = 350.0
 utherm = 20.2678457288
-gamma = 5. / 3.
+gamma = 5.0 / 3.0
 
 start = 0
 stop = 21
 if len(sys.argv) > 1:
-  start = int(sys.argv[1])
+    start = int(sys.argv[1])
 if len(sys.argv) > 2:
-  stop = int(sys.argv[2])
+    stop = int(sys.argv[2])
 
 # Get the analytic solution for the density
 def get_analytic_density(x):
-  return 0.5 * surface_density / scale_height / \
-           np.cosh( (x - x_disc) / scale_height )**2
+    return (
+        0.5 * surface_density / scale_height / np.cosh((x - x_disc) / scale_height) ** 2
+    )
+
 
 # Get the analytic solution for the (isothermal) pressure
 def get_analytic_pressure(x):
-  return (gamma - 1.) * utherm * get_analytic_density(x)
+    return (gamma - 1.0) * utherm * get_analytic_density(x)
+
 
 # Get the data fields to plot from the snapshot file with the given name:
 #  snapshot time, x-coord, density, pressure, velocity norm
 def get_data(name):
-  file = h5py.File(name, "r")
-  coords = np.array(file["/PartType0/Coordinates"])
-  rho = np.array(file["/PartType0/Density"])
-  u = np.array(file["/PartType0/InternalEnergy"])
-  v = np.array(file["/PartType0/Velocities"])
+    file = h5py.File(name, "r")
+    coords = np.array(file["/PartType0/Coordinates"])
+    rho = np.array(file["/PartType0/Density"])
+    u = np.array(file["/PartType0/InternalEnergy"])
+    v = np.array(file["/PartType0/Velocities"])
+
+    P = (gamma - 1.0) * rho * u
 
-  P = (gamma - 1.) * rho * u
+    vtot = np.sqrt(v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2)
 
-  vtot = np.sqrt( v[:,0]**2 + v[:,1]**2 + v[:,2]**2 )
+    return float(file["/Header"].attrs["Time"]), coords[:, 0], rho, P, vtot
 
-  return float(file["/Header"].attrs["Time"]), coords[:,0], rho, P, vtot
 
 # scan the folder for snapshot files and plot all of them (within the requested
 # range)
 for f in sorted(glob.glob("Disc-Patch_*.hdf5")):
-  num = int(f[-8:-5])
-  if num < start or num > stop:
-    continue
-
-  print "processing", f, "..."
-
-  xrange = np.linspace(0., 2. * x_disc, 1000)
-  time, x, rho, P, v = get_data(f)
-
-  fig, ax = pl.subplots(3, 1, sharex = True)
-
-  ax[0].plot(x, rho, "r.")
-  ax[0].plot(xrange, get_analytic_density(xrange), "k-")
-  ax[0].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
-  ax[0].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
-  ax[0].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[0].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[0].set_ylim(0., 1.2 * get_analytic_density(x_disc))
-  ax[0].set_ylabel("density")
-
-  ax[1].plot(x, v, "r.")
-  ax[1].plot(xrange, np.zeros(len(xrange)), "k-")
-  ax[1].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
-  ax[1].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
-  ax[1].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[1].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[1].set_ylim(-0.5, 10.)
-  ax[1].set_ylabel("velocity norm")
-
-  ax[2].plot(x, P, "r.")
-  ax[2].plot(xrange, get_analytic_pressure(xrange), "k-")
-  ax[2].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
-  ax[2].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
-  ax[2].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[2].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
-  ax[2].set_xlim(0., 2. * x_disc)
-  ax[2].set_ylim(0., 1.2 * get_analytic_pressure(x_disc))
-  ax[2].set_xlabel("x")
-  ax[2].set_ylabel("pressure")
-
-  pl.suptitle("t = {0:.2f}".format(time))
-
-  pl.savefig("{name}.png".format(name = f[:-5]))
-  pl.close()
+    num = int(f[-8:-5])
+    if num < start or num > stop:
+        continue
+
+    print("processing", f, "...")
+
+    xrange = np.linspace(0.0, 2.0 * x_disc, 1000)
+    time, x, rho, P, v = get_data(f)
+
+    fig, ax = pl.subplots(3, 1, sharex=True)
+
+    ax[0].plot(x, rho, "r.")
+    ax[0].plot(xrange, get_analytic_density(xrange), "k-")
+    ax[0].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+    ax[0].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+    ax[0].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[0].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[0].set_ylim(0.0, 1.2 * get_analytic_density(x_disc))
+    ax[0].set_ylabel("density")
+
+    ax[1].plot(x, v, "r.")
+    ax[1].plot(xrange, np.zeros(len(xrange)), "k-")
+    ax[1].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+    ax[1].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+    ax[1].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[1].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[1].set_ylim(-0.5, 10.0)
+    ax[1].set_ylabel("velocity norm")
+
+    ax[2].plot(x, P, "r.")
+    ax[2].plot(xrange, get_analytic_pressure(xrange), "k-")
+    ax[2].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+    ax[2].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+    ax[2].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[2].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+    ax[2].set_xlim(0.0, 2.0 * x_disc)
+    ax[2].set_ylim(0.0, 1.2 * get_analytic_pressure(x_disc))
+    ax[2].set_xlabel("x")
+    ax[2].set_ylabel("pressure")
+
+    pl.suptitle("t = {0:.2f}".format(time))
+
+    pl.savefig("{name}.png".format(name=f[:-5]))
+    pl.close()

examples/GravityTests/DiscPatch/HydroStatic_1D/run.sh

@@ -4,10 +4,10 @@
 if [ ! -e Disc-Patch.hdf5 ]
 then
     echo "Generating initial conditions for the disc patch example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../../swift --external-gravity --hydro --threads=4 disc-patch-icc.yml 2>&1 | tee output.log
+../../../../swift --external-gravity --hydro --threads=4 disc-patch-icc.yml 2>&1 | tee output.log
 
-python plotSolution.py
+python3 plotSolution.py

examples/GravityTests/ExternalPointMass/energy_plot.py

 import matplotlib
+
 matplotlib.use("Agg")
 from pylab import *
 import h5py
 
 # Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 12,
-'legend.fontsize': 12,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (3.15,3.15),
-'figure.subplot.left'    : 0.145,
-'figure.subplot.right'   : 0.99,
-'figure.subplot.bottom'  : 0.11,
-'figure.subplot.top'     : 0.99,
-'figure.subplot.wspace'  : 0.15,
-'figure.subplot.hspace'  : 0.12,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
+params = {
+    "axes.labelsize": 10,
+    "axes.titlesize": 10,
+    "font.size": 12,
+    "legend.fontsize": 12,
+    "xtick.labelsize": 10,
+    "ytick.labelsize": 10,
+    "text.usetex": True,
+    "figure.figsize": (3.15, 3.15),
+    "figure.subplot.left": 0.145,
+    "figure.subplot.right": 0.99,
+    "figure.subplot.bottom": 0.11,
+    "figure.subplot.top": 0.99,
+    "figure.subplot.wspace": 0.15,
+    "figure.subplot.hspace": 0.12,
+    "lines.markersize": 6,
+    "lines.linewidth": 3.0,
 }
 rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
 
 
 import numpy as np
@@ -35,7 +35,7 @@ stats_filename = "./energy.txt"
 
 # First snapshot
 snap_filename = "pointMass_0000.hdf5"
-f = h5.File(snap_filename,'r')
+f = h5.File(snap_filename, "r")
 
 # Read the units parameters from the snapshot
 units = f["InternalCodeUnits"]
@@ -43,7 +43,7 @@ unit_mass = units.attrs["Unit mass in cgs (U_M)"]
 unit_length = units.attrs["Unit length in cgs (U_L)"]
 unit_time = units.attrs["Unit time in cgs (U_t)"]
 
-G = 6.67408e-8 * unit_mass * unit_time**2 / unit_length**3
+G = 6.67408e-8 * unit_mass * unit_time ** 2 / unit_length ** 3
 
 # Read the header
 header = f["Header"]
@@ -52,14 +52,14 @@ box_size = float(header.attrs["BoxSize"][0])
 # Read the properties of the potential
 parameters = f["Parameters"]
 mass = float(parameters.attrs["PointMassPotential:mass"])
-centre = [box_size/2, box_size/2, box_size/2]
+centre = [box_size / 2, box_size / 2, box_size / 2]
 f.close()
 
 # Read the statistics summary
 file_energy = np.loadtxt("energy.txt")
-time_stats = file_energy[:,0]
-E_kin_stats = file_energy[:,3]
-E_pot_stats = file_energy[:,5]
+time_stats = file_energy[:, 0]
+E_kin_stats = file_energy[:, 3]
+E_pot_stats = file_energy[:, 5]
 E_tot_stats = E_kin_stats + E_pot_stats
 
 # Read the snapshots
@@ -71,29 +71,33 @@ Lz_snap = np.zeros(402)
 
 # Read all the particles from the snapshots
 for i in range(402):
-    snap_filename = "pointMass_%0.4d.hdf5"%i
-    f = h5.File(snap_filename,'r')
-
-    pos_x = f["PartType1/Coordinates"][:,0]
-    pos_y = f["PartType1/Coordinates"][:,1]
-    pos_z = f["PartType1/Coordinates"][:,2]
-    vel_x = f["PartType1/Velocities"][:,0]
-    vel_y = f["PartType1/Velocities"][:,1]
-    vel_z = f["PartType1/Velocities"][:,2]
+    snap_filename = "pointMass_%0.4d.hdf5" % i
+    f = h5.File(snap_filename, "r")
+
+    pos_x = f["PartType1/Coordinates"][:, 0]
+    pos_y = f["PartType1/Coordinates"][:, 1]
+    pos_z = f["PartType1/Coordinates"][:, 2]
+    vel_x = f["PartType1/Velocities"][:, 0]
+    vel_y = f["PartType1/Velocities"][:, 1]
+    vel_z = f["PartType1/Velocities"][:, 2]
     m = f["/PartType1/Masses"][:]
-    
-    r = np.sqrt((pos_x[:] - centre[0])**2 + (pos_y[:] - centre[1])**2 + (pos_z[:] - centre[2])**2)
+
+    r = np.sqrt(
+        (pos_x[:] - centre[0]) ** 2
+        + (pos_y[:] - centre[1]) ** 2
+        + (pos_z[:] - centre[2]) ** 2
+    )
     Lz = (pos_x[:] - centre[0]) * vel_y[:] - (pos_y[:] - centre[1]) * vel_x[:]
 
     time_snap[i] = f["Header"].attrs["Time"]
-    E_kin_snap[i] = np.sum(0.5 * m * (vel_x[:]**2 + vel_y[:]**2 + vel_z[:]**2))
+    E_kin_snap[i] = np.sum(0.5 * m * (vel_x[:] ** 2 + vel_y[:] ** 2 + vel_z[:] ** 2))
     E_pot_snap[i] = np.sum(-mass * m * G / r)
     E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
     Lz_snap[i] = np.sum(Lz)
 
 print("Starting energy:", E_kin_stats[0], E_pot_stats[0], E_tot_stats[0])
 print("Ending   energy:", E_kin_stats[-1], E_pot_stats[-1], E_tot_stats[-1])
-    
+
 # Plot energy evolution
 figure()
 plot(time_stats, E_kin_stats, "r-", lw=0.5, label="Kinetic energy")
@@ -106,7 +110,7 @@ plot(time_snap[::10], E_tot_snap[::10], "kD", lw=0.5, ms=2)
 
 legend(loc="center right", fontsize=8, frameon=False, handlelength=3, ncol=1)
 xlabel("${\\rm{Time}}$", labelpad=0)
-ylabel("${\\rm{Energy}}$",labelpad=0)
+ylabel("${\\rm{Energy}}$", labelpad=0)
 xlim(0, 8)
 
 savefig("energy.png", dpi=200)
@@ -116,9 +120,7 @@ figure()
 plot(time_snap, Lz_snap, "k-", lw=0.5, ms=2)
 
 xlabel("${\\rm{Time}}$", labelpad=0)
-ylabel("${\\rm{Angular~momentum}}$",labelpad=0)
+ylabel("${\\rm{Angular~momentum}}$", labelpad=0)
 xlim(0, 8)
 
 savefig("angular_momentum.png", dpi=200)
-
-

examples/GravityTests/ExternalPointMass/makeIC.py

 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
- #                    Tom Theuns (tom.theuns@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/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
+#                    Tom Theuns (tom.theuns@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 h5py
 import sys
@@ -27,50 +27,55 @@ import random
 # Generates a random distriution of particles, for motion in an external potential centred at (0,0,0)
 
 # physical constants in cgs
-NEWTON_GRAVITY_CGS  = 6.67408e-8
-SOLAR_MASS_IN_CGS   = 1.98848e33
-PARSEC_IN_CGS       = 3.08567758e18
+NEWTON_GRAVITY_CGS = 6.67408e-8
+SOLAR_MASS_IN_CGS = 1.98848e33
+PARSEC_IN_CGS = 3.08567758e18
 
 # choice of units
-const_unit_length_in_cgs   =   (1000*PARSEC_IN_CGS)
-const_unit_mass_in_cgs     =   (SOLAR_MASS_IN_CGS)
-const_unit_velocity_in_cgs =   (1e5)
+const_unit_length_in_cgs = 1000 * PARSEC_IN_CGS
+const_unit_mass_in_cgs = SOLAR_MASS_IN_CGS
+const_unit_velocity_in_cgs = 1e5
 
-print("UnitMass_in_cgs:     ", const_unit_mass_in_cgs) 
+print("UnitMass_in_cgs:     ", const_unit_mass_in_cgs)
 print("UnitLength_in_cgs:   ", const_unit_length_in_cgs)
 print("UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs)
 print("UnitTime_in_cgs:     ", const_unit_length_in_cgs / const_unit_velocity_in_cgs)
 
 # derived units
-const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
-const_G                = ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
-print('---------------------')
-print('G in internal units: ', const_G)
+const_unit_time_in_cgs = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+const_G = (
+    NEWTON_GRAVITY_CGS
+    * const_unit_mass_in_cgs
+    * const_unit_time_in_cgs
+    * const_unit_time_in_cgs
+    / (const_unit_length_in_cgs * const_unit_length_in_cgs * const_unit_length_in_cgs)
+)
+print("---------------------")
+print("G in internal units: ", const_G)
 
 
 # Parameters
-periodic   = 1            # 1 For periodic box
-boxSize    = 100.         # 
-max_radius = boxSize / 4. # maximum radius of particles
-Mass       = 1e10         
+periodic = 1  # 1 For periodic box
+boxSize = 100.0  #
+max_radius = boxSize / 4.0  # maximum radius of particles
+Mass = 1e10
 print("Mass at the centre:  ", Mass)
 
 numPart = int(sys.argv[1])  # Number of particles
-mass    = 1.
+mass = 1.0
 
-fileName = "PointMass.hdf5" 
+fileName = "PointMass.hdf5"
 
 
+# --------------------------------------------------
 
-#--------------------------------------------------
-
-#File
-file = h5py.File(fileName, 'w')
+# File
+file = h5py.File(fileName, "w")
 
 # Header
 grp = file.create_group("/Header")
 grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [0, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] = [0, numPart, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [0, numPart, 0, 0, 0, 0]
 grp.attrs["Time"] = 0.0
@@ -80,52 +85,52 @@ grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["Dimension"] = 3
 
 
-#Units
+# Units
 grp = file.create_group("/Units")
 grp.attrs["Unit length in cgs (U_L)"] = 3.0856776e21
 grp.attrs["Unit mass in cgs (U_M)"] = 1.98848e33
 grp.attrs["Unit time in cgs (U_t)"] = 3.0856776e16
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
-#Particle group
+# Particle group
 grp1 = file.create_group("/PartType1")
 
-#generate particle positions
-radius = max_radius * (numpy.random.rand(numPart))**(1./3.)
-print('---------------------')
-print('Radius: minimum = ',min(radius))
-print('Radius: maximum = ',max(radius))
+# generate particle positions
+radius = max_radius * (numpy.random.rand(numPart)) ** (1.0 / 3.0)
+print("---------------------")
+print("Radius: minimum = ", min(radius))
+print("Radius: maximum = ", max(radius))
 radius = numpy.sort(radius)
-r      = numpy.zeros((numPart, 3))
-r[:,0] = radius
-
-#generate particle velocities
-speed  = numpy.sqrt(const_G * Mass / radius)
-omega  = speed / radius
-period = 2.*math.pi/omega
-print('---------------------')
-print('Period: minimum = ',min(period))
-print('Period: maximum = ',max(period))
-
-v      = numpy.zeros((numPart, 3))
-v[:,0] = -omega * r[:,1]
-v[:,1] =  omega * r[:,0]
-
-ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+r = numpy.zeros((numPart, 3))
+r[:, 0] = radius
+
+# generate particle velocities
+speed = numpy.sqrt(const_G * Mass / radius)
+omega = speed / radius
+period = 2.0 * math.pi / omega
+print("---------------------")
+print("Period: minimum = ", min(period))
+print("Period: maximum = ", max(period))
+
+v = numpy.zeros((numPart, 3))
+v[:, 0] = -omega * r[:, 1]
+v[:, 1] = omega * r[:, 0]
+
+ds = grp1.create_dataset("Velocities", (numPart, 3), "f")
 ds[()] = v
 v = numpy.zeros(1)
 
-m = numpy.full((numPart, ), mass)
-ds = grp1.create_dataset('Masses', (numPart,), 'f')
+m = numpy.full((numPart,), mass)
+ds = grp1.create_dataset("Masses", (numPart,), "f")
 ds[()] = m
 m = numpy.zeros(1)
 
 ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
-ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
+ds = grp1.create_dataset("ParticleIDs", (numPart,), "L")
 ds[()] = ids
 
-ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ds = grp1.create_dataset("Coordinates", (numPart, 3), "d")
 ds[()] = r
 
 file.close()

examples/GravityTests/ExternalPointMass/run.sh

@@ -4,10 +4,10 @@
 if [ ! -e PointMass.hdf5 ]
 then
     echo "Generating initial conditions for the point mass potential box example..."
-    python makeIC.py 10000
+    python3 makeIC.py 10000
 fi
 
 rm -rf pointMass_*.hdf5
-../../swift --external-gravity --threads=1 externalPointMass.yml 2>&1 | tee output.log
+../../../swift --external-gravity --threads=1 externalPointMass.yml 2>&1 | tee output.log
 
-python energy_plot.py
+python3 energy_plot.py

examples/GravityTests/Gravity_glass/makeIC.py

 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+#                    Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# 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 h5py
 import sys
@@ -26,25 +26,25 @@ import numpy as np
 # with a density of 1
 
 # Parameters
-periodic= 1           # 1 For periodic box
-boxSize = 1.
-rho = 1.
+periodic = 1  # 1 For periodic box
+boxSize = 1.0
+rho = 1.0
 L = int(sys.argv[1])  # Number of particles along one axis
 fileName = "uniform_DM_box.hdf5"
 
-#---------------------------------------------------
-numPart = L**3
-mass = boxSize**3 * rho / numPart
+# ---------------------------------------------------
+numPart = L ** 3
+mass = boxSize ** 3 * rho / numPart
 
-#--------------------------------------------------
+# --------------------------------------------------
 
-#File
-file = h5py.File(fileName, 'w')
+# File
+file = h5py.File(fileName, "w")
 
 # Header
 grp = file.create_group("/Header")
 grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [0, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] = [0, numPart, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [0, numPart, 0, 0, 0, 0]
 grp.attrs["Time"] = 0.0
@@ -53,28 +53,28 @@ grp.attrs["MassTable"] = [0.0, mass, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = 0
 grp.attrs["Dimension"] = 3
 
-#Units
+# Units
 grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.
-grp.attrs["Unit mass in cgs (U_M)"] = 1.
-grp.attrs["Unit time in cgs (U_t)"] = 1.
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+grp.attrs["Unit length in cgs (U_L)"] = 1.0
+grp.attrs["Unit mass in cgs (U_M)"] = 1.0
+grp.attrs["Unit time in cgs (U_t)"] = 1.0
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
 
-#Particle group
+# Particle group
 grp = file.create_group("/PartType1")
 
-v  = np.zeros((numPart, 3))
-ds = grp.create_dataset('Velocities', (numPart, 3), 'f', data=v)
+v = np.zeros((numPart, 3))
+ds = grp.create_dataset("Velocities", (numPart, 3), "f", data=v)
 
 m = np.full((numPart, 1), mass)
-ds = grp.create_dataset('Masses', (numPart,1), 'f', data=m)
+ds = grp.create_dataset("Masses", (numPart, 1), "f", data=m)
 
-ids = np.linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
-ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ids = np.linspace(0, numPart, numPart, endpoint=False).reshape((numPart, 1))
+ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
 ds[()] = ids + 1
 coords = np.random.rand(numPart, 3) * boxSize
-ds = grp.create_dataset('Coordinates', (numPart, 3), 'd', data=coords)
+ds = grp.create_dataset("Coordinates", (numPart, 3), "d", data=coords)
 
 file.close()

examples/GravityTests/Hernquist_circularorbit/makeIC.py

@@ -16,9 +16,6 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 ################################################################################
-from galpy.potential import NFWPotential
-from galpy.orbit import Orbit
-from galpy.util import bovy_conversion
 import numpy as np
 import matplotlib.pyplot as plt
 from astropy import units
@@ -27,7 +24,7 @@ import h5py as h5
 C = 8.0
 M_200 = 2.0
 N_PARTICLES = 3
-print("Initial conditions written to 'test_nfw.hdf5'")
+print("Initial conditions written to 'circularorbitshernquist.hdf5'")
 
 pos = np.zeros((3, 3))
 pos[0, 2] = 50.0

examples/GravityTests/Hernquist_circularorbit/plotprog.py

@@ -22,7 +22,7 @@ import h5py
 import matplotlib.pyplot as plt
 from scipy.integrate import odeint
 
-t = np.linspace(0, 40, 1e5)
+t = np.linspace(0, 40, 100000)
 y0 = [0, 10]
 a = 30.0
 G = 4.300927e-06

examples/GravityTests/Hernquist_circularorbit/run.sh

@@ -6,18 +6,18 @@ then
     if command -v python3 &>/dev/null; then
         python3 makeIC.py
     else 
-        python makeIC.py
+        python3 makeIC.py
     fi
 
 fi
 
 # self gravity G, external potential g, hydro s, threads t and high verbosity v
-../../swift --external-gravity --threads=6 hernquistcirc.yml 2>&1 | tee output.log
+../../../swift --external-gravity --threads=6 hernquistcirc.yml 2>&1 | tee output.log
 
 
 echo "Save plots of the circular orbits"
 if command -v python3 &>/dev/null; then
     python3 plotprog.py
 else 
-    python plotprog.py
+    python3 plotprog.py
 fi

examples/GravityTests/Hernquist_radialinfall/makeIC.py

@@ -25,8 +25,8 @@ import random
 import numpy as np
 
 # Generates N particles in a spherical distribution centred on [0,0,0], to be moved in an isothermal potential
-# usage: python makeIC.py 1000 0 : generate 1000 particles on circular orbits
-#        python makeIC.py 1000 1 : generate 1000 particles with Lz/L uniform in [0,1]
+# usage: python3 makeIC.py 1000 0 : generate 1000 particles on circular orbits
+#        python3 makeIC.py 1000 1 : generate 1000 particles with Lz/L uniform in [0,1]
 # all particles move in the xy plane, and start at y=0
 
 # physical constants in cgs

examples/GravityTests/Hernquist_radialinfall/run.sh

@@ -7,12 +7,12 @@ then
     if command -v python3 &>/dev/null; then
         python3 makeIC.py 
     else 
-        python makeIC.py
+        python3 makeIC.py
     fi
 fi
 
 rm -rf hernquist_*.hdf5
-../../swift --external-gravity --threads=1 hernquist.yml 2>&1 | tee output.log
+../../../swift --external-gravity --threads=1 hernquist.yml 2>&1 | tee output.log
 
 
 
@@ -20,5 +20,5 @@ echo "Make plots of the radially free falling particles"
 if command -v python3 &>/dev/null; then
     python3 plotprog.py 
 else 
-    python plotprog.py
+    python3 plotprog.py
 fi

examples/GravityTests/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/>.
- # 
- ##############################################################################
+# 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
+
 matplotlib.use("Agg")
 from pylab import *
 import sys
 
-#for the plotting
-max_r = float(sys.argv[1]) #in units of the virial radius 
+# for the plotting
+max_r = float(sys.argv[1])  # in units of the virial radius
 n_radial_bins = int(sys.argv[2])
 n_snaps = int(sys.argv[3])
 
-#some constants
-OMEGA = 0.3 # Cosmological matter fraction at z = 0
+# some constants
+OMEGA = 0.3  # Cosmological matter fraction at z = 0
 PARSEC_IN_CGS = 3.0856776e18
 KM_PER_SEC_IN_CGS = 1.0e5
 CONST_G_CGS = 6.672e-8
 CONST_m_H_CGS = 1.67e-24
-h = 0.67777 # hubble parameter
-gamma = 5./3.
+h = 0.67777  # hubble parameter
+gamma = 5.0 / 3.0
 eta = 1.2349
-H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
 
-#read some header/parameter information from the first snapshot
+# read some header/parameter information from the first snapshot
 
 filename = "Hydrostatic_0000.hdf5"
-f = h5.File(filename,'r')
+f = h5.File(filename, "r")
 params = f["Parameters"]
 unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
@@ -51,80 +52,86 @@ 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"])
+# 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"])
 
-#calculate r_vir and M_vir from v_c
-r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
-M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+# calculate r_vir and M_vir from v_c
+r_vir_cgs = v_c_cgs / (10.0 * H_0_cgs * np.sqrt(OMEGA))
+M_vir_cgs = r_vir_cgs * v_c_cgs ** 2 / CONST_G_CGS
 
 for i in range(n_snaps):
 
-    filename = "Hydrostatic_%04d.hdf5" %i
-    f = h5.File(filename,'r')
+    filename = "Hydrostatic_%04d.hdf5" % i
+    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
-    coords[:,0] -= box_centre[0]/2.
-    coords[:,1] -= box_centre[1]/2.
-    coords[:,2] -= box_centre[2]/2.
-    radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2)
-    radius_cgs = radius*unit_length_cgs
+    coords[:, 0] -= box_centre[0] / 2.0
+    coords[:, 1] -= box_centre[1] / 2.0
+    coords[:, 2] -= box_centre[2] / 2.0
+    radius = np.sqrt(coords[:, 0] ** 2 + coords[:, 1] ** 2 + coords[:, 2] ** 2)
+    radius_cgs = radius * unit_length_cgs
     radius_over_virial_radius = radius_cgs / r_vir_cgs
 
     r = radius_over_virial_radius
 
-    bin_edges = np.linspace(0.,max_r,n_radial_bins+1)
+    bin_edges = np.linspace(0.0, max_r, n_radial_bins + 1)
     bin_width = bin_edges[1] - bin_edges[0]
-    hist = np.histogram(r,bins = bin_edges)[0] # number of particles in each bin
-
+    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 
+    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 = 4.*np.pi * radial_bin_mids**2 * bin_width
+    radial_bin_mids = np.linspace(
+        bin_width / 2.0, max_r - bin_width / 2.0, n_radial_bins
+    )
 
+    # volume in each radial bin
+    volume = 4.0 * 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
 
-    t = np.linspace(10./n_radial_bins,10.0,1000)
-    rho_analytic = t**(-2)/(4.*np.pi)
+    t = np.linspace(10.0 / n_radial_bins, 10.0, 1000)
+    rho_analytic = t ** (-2) / (4.0 * np.pi)
 
-
-    #initial analytic density profile    
-    if (i == 0):
+    # 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)
+        rho_analytic_init = rho_0 * (radial_bin_mids / r_0) ** (-2)
 
     figure()
-    plot(radial_bin_mids,density/rho_analytic_init,'ko',label = "Average density of shell")
-    #plot(t,rho_analytic,label = "Initial analytic density profile")
+    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
-    savefig(plot_filename,format = "png")
+    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
+    savefig(plot_filename, format="png")
     close()
-

examples/GravityTests/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/>.
- # 
- ##############################################################################
+# 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
+
 matplotlib.use("Agg")
 from pylab import *
 import sys
 
-def do_binning(x,y,x_bin_edges):
 
-    #x and y are arrays, where y = f(x)
-    #returns number of elements of x in each bin, and the total of the y elements corresponding to those x values
+def do_binning(x, y, x_bin_edges):
+
+    # x and y are arrays, where y = f(x)
+    # returns number of elements of x in each bin, and the total of the y elements corresponding to those x values
 
     n_bins = x_bin_edges.size - 1
     count = np.zeros(n_bins)
     y_totals = np.zeros(n_bins)
-    
+
     for i in range(n_bins):
-        ind = np.intersect1d(np.where(x > bin_edges[i])[0],np.where(x <= bin_edges[i+1])[0])
+        ind = np.intersect1d(
+            np.where(x > bin_edges[i])[0], np.where(x <= bin_edges[i + 1])[0]
+        )
         count[i] = ind.size
         binned_y = y[ind]
         y_totals[i] = np.sum(binned_y)
 
-    return(count,y_totals)
+    return (count, y_totals)
 
 
-#for the plotting
+# for the plotting
 max_r = float(sys.argv[1])
 n_radial_bins = int(sys.argv[2])
 n_snaps = int(sys.argv[3])
 
-#some constants
-OMEGA = 0.3 # Cosmological matter fraction at z = 0
+# some constants
+OMEGA = 0.3  # Cosmological matter fraction at z = 0
 PARSEC_IN_CGS = 3.0856776e18
 KM_PER_SEC_IN_CGS = 1.0e5
 CONST_G_CGS = 6.672e-8
 CONST_m_H_CGS = 1.67e-24
-h = 0.67777 # hubble parameter
-gamma = 5./3.
+h = 0.67777  # hubble parameter
+gamma = 5.0 / 3.0
 eta = 1.2349
-H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
 
-#read some header/parameter information from the first snapshot
+# read some header/parameter information from the first snapshot
 
 filename = "Hydrostatic_0000.hdf5"
-f = h5.File(filename,'r')
+f = h5.File(filename, "r")
 params = f["Parameters"]
 unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
@@ -73,54 +77,55 @@ header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
 box_centre = np.array(header.attrs["BoxSize"])
 
-#calculate r_vir and M_vir from v_c
-r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
-M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+# calculate r_vir and M_vir from v_c
+r_vir_cgs = v_c_cgs / (10.0 * H_0_cgs * np.sqrt(OMEGA))
+M_vir_cgs = r_vir_cgs * v_c_cgs ** 2 / CONST_G_CGS
 
 for i in range(n_snaps):
 
-    filename = "Hydrostatic_%04d.hdf5" %i
-    f = h5.File(filename,'r')
+    filename = "Hydrostatic_%04d.hdf5" % i
+    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
-    coords[:,0] -= box_centre[0]/2.
-    coords[:,1] -= box_centre[1]/2.
-    coords[:,2] -= box_centre[2]/2.
-    radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2)
-    radius_cgs = radius*unit_length_cgs
+    coords[:, 0] -= box_centre[0] / 2.0
+    coords[:, 1] -= box_centre[1] / 2.0
+    coords[:, 2] -= box_centre[2] / 2.0
+    radius = np.sqrt(coords[:, 0] ** 2 + coords[:, 1] ** 2 + coords[:, 2] ** 2)
+    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
-    u /= v_c**2/(2. * (gamma - 1.))
+    # make dimensionless
+    u /= v_c ** 2 / (2.0 * (gamma - 1.0))
     r = radius_over_virial_radius
 
-    bin_edges = np.linspace(0,max_r,n_radial_bins + 1)
-    (hist,u_totals) = do_binning(r,u,bin_edges)
-    
+    bin_edges = np.linspace(0, max_r, n_radial_bins + 1)
+    (hist, u_totals) = do_binning(r, u, bin_edges)
+
     bin_widths = bin_edges[1] - bin_edges[0]
-    radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins) 
+    radial_bin_mids = np.linspace(
+        bin_widths / 2.0, max_r - bin_widths / 2.0, n_radial_bins
+    )
     binned_u = u_totals / hist
 
     figure()
-    plot(radial_bin_mids,binned_u,'ko',label = "Numerical solution")
-    legend(loc = "lower right")
+    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
-    savefig(plot_filename,format = "png")
+    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
+    savefig(plot_filename, format="png")
     close()
-
-
-        
-