diff --git a/examples/Disk-Patch/GravityOnly/README b/examples/Disk-Patch/GravityOnly/README
new file mode 100644
index 0000000000000000000000000000000000000000..5bf2638fc5ed48ebe248773223dde888af0c3bc8
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/README
@@ -0,0 +1,10 @@
+Setup for a potential of a patch disk, see Creasey, Theuns &
+Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
+
+The density is given by
+rho(z) = (Sigma/2b) / cosh^2(z/b)
+where Sigma is the surface density, and b the scale height.
+
+The corresponding force is
+dphi/dz = 2 pi G Sigma tanh(z/b),
+which satifies d^2phi/dz^2 = 4 pi G rho.
diff --git a/examples/Disk-Patch/GravityOnly/disk-patch.yml b/examples/Disk-Patch/GravityOnly/disk-patch.yml
new file mode 100644
index 0000000000000000000000000000000000000000..78b42e78356f83e80eee8e7f5f91ad7dcf90c37f
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/disk-patch.yml
@@ -0,0 +1,43 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.9885e33 # Grams
+ UnitLength_in_cgs: 3.0856776e18 # Centimeters
+ UnitVelocity_in_cgs: 1e5 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 480. # The end time of the simulation (in internal units).
+ dt_min: 1e-3 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1.0 # Time between statistics output
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: Disk-Patch # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 8. # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 1. # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
+ max_smoothing_length: 40. # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: Disk-Patch.hdf5 # The file to read
+
+# External potential parameters
+Disk-PatchPotential:
+ surface_density: 10.
+ scale_height: 100.
+ z_disk: 300.
+ timestep_mult: 0.03
diff --git a/examples/Disk-Patch/GravityOnly/makeIC.py b/examples/Disk-Patch/GravityOnly/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..135d1d4e34e0ec4c32d75a2c49d6e55f46db4e1c
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/makeIC.py
@@ -0,0 +1,161 @@
+###############################################################################
+ # 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
+import numpy
+import math
+import random
+
+# Generates N particles in a box of [0:BoxSize,0:BoxSize,-2scale_height:2scale_height]
+# see Creasey, Theuns & Bower, 2013, for the equations:
+# disk parameters are: surface density sigma
+# scale height b
+# density: rho(z) = (sigma/2b) sech^2(z/b)
+# isothermal velocity dispersion = <v_z^2? = b pi G sigma
+# grad potential = 2 pi G sigma tanh(z/b)
+# potential = ln(cosh(z/b)) + const
+# Dynamical time = sqrt(b / (G sigma))
+# to obtain the 1/ch^2(z/b) profile from a uniform profile (a glass, say, or a uniform random variable), note that, when integrating in z
+# \int 0^z dz/ch^2(z) = tanh(z)-tanh(0) = \int_0^x dx = x (where the last integral refers to a uniform density distribution), so that z = atanh(x)
+# usage: python makeIC.py 1000
+
+# physical constants in cgs
+NEWTON_GRAVITY_CGS = 6.672e-8
+SOLAR_MASS_IN_CGS = 1.9885e33
+PARSEC_IN_CGS = 3.0856776e18
+PROTON_MASS_IN_CGS = 1.6726231e24
+YEAR_IN_CGS = 3.154e+7
+
+# 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)
+
+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
+
+
+# parameters of potential
+surface_density = 10.
+scale_height = 100.
+
+# 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 'G=', const_G
+v_disp = numpy.sqrt(scale_height * math.pi * const_G * surface_density)
+t_dyn = numpy.sqrt(scale_height / (const_G * surface_density))
+print 'dynamical time = ',t_dyn
+print ' velocity dispersion = ',v_disp
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 600. #
+Radius = 100. # maximum radius of particles [kpc]
+G = const_G
+
+N = int(sys.argv[1]) # Number of particles
+
+# these are not used but necessary for I/O
+rho = 2. # Density
+P = 1. # Pressure
+gamma = 5./3. # Gas adiabatic index
+fileName = "Disk-Patch.hdf5"
+
+
+#---------------------------------------------------
+numPart = N
+mass = 1
+internalEnergy = P / ((gamma - 1.)*rho)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs
+grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+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
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+# set seed for random number
+numpy.random.seed(1234)
+
+#Particle group
+#grp0 = file.create_group("/PartType0")
+grp1 = file.create_group("/PartType1")
+
+#generate particle positions
+r = numpy.zeros((numPart, 3))
+r[:,0] = numpy.random.rand(N) * boxSize
+r[:,1] = numpy.random.rand(N) * boxSize
+z = scale_height * numpy.arctanh(numpy.random.rand(2*N))
+gd = z < boxSize / 2
+r[:,2] = z[gd][0:N]
+random = numpy.random.rand(N) > 0.5
+r[random,2] *= -1
+r[:,2] += 0.5 * boxSize
+
+#generate particle velocities
+v = numpy.zeros((numPart, 3))
+v = numpy.zeros(1)
+#v[:,2] =
+
+
+ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+
+
+m = numpy.ones((numPart, ), dtype=numpy.float32) * mass
+ds = grp1.create_dataset('Masses', (numPart,), 'f')
+ds[()] = m
+m = numpy.zeros(1)
+
+
+ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
+ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
+ds[()] = ids
+
+ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = r
+
+
+file.close()
diff --git a/examples/Disk-Patch/GravityOnly/run.sh b/examples/Disk-Patch/GravityOnly/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a123ad24d7ca34105c22f5f31e75c688c681288f
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e Isothermal.hdf5 ]
+then
+ echo "Generating initial conditions for the disk-patch example..."
+ python makeIC.py 1000
+fi
+
+../../swift -g -t 2 disk-patch.yml
diff --git a/examples/Disk-Patch/GravityOnly/test.pro b/examples/Disk-Patch/GravityOnly/test.pro
new file mode 100644
index 0000000000000000000000000000000000000000..4bd0d001975d80b6729cf2ef7b95f81da5bc4fe8
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/test.pro
@@ -0,0 +1,158 @@
+;
+; 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 = 'Disk-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
+
+
diff --git a/examples/Disk-Patch/HydroStatic/README b/examples/Disk-Patch/HydroStatic/README
new file mode 100644
index 0000000000000000000000000000000000000000..1bbb7e8a34f261e0feb1cab60d9d31519817b0b2
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/README
@@ -0,0 +1,19 @@
+Generates and evolves a disk-patch, where gas is in hydrostatic
+equilibrium with an imposed external gravitational force, using the
+equations from Creasey, Theuns & Bower, 2013, MNRAS, Volume 429,
+Issue 3, p.1922-1948.
+
+To generate ICs ready for a scientific run:
+
+1) Recover a uniform glass file by running 'getGlass.sh'.
+
+2) Generate pre-ICs by running the 'makeIC.py' script.
+
+3) Run SWIFT with an isothermal EoS, no cooling nor feedback, and the
+disk-patch potential switched on and using the parameters from
+'disk-patch-icc.yml'
+
+4) The ICs are then ready to be run for a science problem.
+
+When running SWIFT with the parameters from 'disk-patch.yml' and an
+ideal gas EoS on these ICs the disk should stay in equilibrium.
diff --git a/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml b/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml
new file mode 100644
index 0000000000000000000000000000000000000000..ebf04675852a7663119ed1ecfd33a05da6b7bb15
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml
@@ -0,0 +1,44 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.9885e33 # Grams
+ UnitLength_in_cgs: 3.0856776e18 # Centimeters
+ UnitVelocity_in_cgs: 1e5 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0 # The starting time of the simulation (in internal units).
+ time_end: 968. # The end time of the simulation (in internal units).
+ dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1. # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1 # Time between statistics output
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: Disk-Patch # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 12. # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
+ max_smoothing_length: 70. # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: Disk-Patch.hdf5 # The file to read
+
+# External potential parameters
+Disk-PatchPotential:
+ surface_density: 10.
+ scale_height: 100.
+ z_disk: 200.
+ timestep_mult: 0.03
+ growth_time: 5.
diff --git a/examples/Disk-Patch/HydroStatic/disk-patch.yml b/examples/Disk-Patch/HydroStatic/disk-patch.yml
new file mode 100644
index 0000000000000000000000000000000000000000..55df81a08d16c6a4f39aa5e9e9205101dedaa3a9
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/disk-patch.yml
@@ -0,0 +1,43 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.9885e33 # Grams
+ UnitLength_in_cgs: 3.0856776e18 # Centimeters
+ UnitVelocity_in_cgs: 1e5 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 968 # The starting time of the simulation (in internal units).
+ time_end: 12000. # The end time of the simulation (in internal units).
+ dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1. # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1 # Time between statistics output
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: Disk-Patch-dynamic # Common part of the name of output files
+ time_first: 968. # Time of the first output (in internal units)
+ delta_time: 24. # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
+ max_smoothing_length: 70. # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: Disk-Patch-dynamic.hdf5 # The file to read
+
+# External potential parameters
+Disk-PatchPotential:
+ surface_density: 10.
+ scale_height: 100.
+ z_disk: 200.
+ timestep_mult: 0.03
diff --git a/examples/Disk-Patch/HydroStatic/dynamic.pro b/examples/Disk-Patch/HydroStatic/dynamic.pro
new file mode 100644
index 0000000000000000000000000000000000000000..c02c65fe418e84cdd62978dbddcf5a641fa4c156
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/dynamic.pro
@@ -0,0 +1,129 @@
+;
+; 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 = 'Disk-Patch-dynamic_'
+
+; 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.
+z_disk = 200.;
+gamma = 5./3.
+
+; derived units
+constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
+pcentre = [0.,0.,z_disk] * pc / uL
+utherm = !pi * constG * surface_density * scale_height / (gamma-1.)
+soundspeed = sqrt(gamma * (gamma-1.) * utherm)
+t_dyn = sqrt(scale_height / (constG * surface_density))
+
+;
+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,'PartType0/ParticleIDs')
+nfollow = n_elements(id)
+
+
+; compute anlytic profile
+nbins = 100
+zbins = findgen(nbins)/float(nbins-1) * 2 * scale_height
+rbins = (surface_density/(2.*scale_height)) / cosh(abs(zbins)/scale_height)^2
+
+
+; plot analytic profile
+wset,0
+plot,[0],[0],xr=[0,2*scale_height],yr=[0,max(rbins)],/nodata,xtitle='|z|',ytitle=textoidl('\rho')
+oplot,zbins,rbins,color=blue
+
+ifile = 0
+nskip = nfiles - 1
+isave = 0
+nplot = 8192 ; randomly plot particles
+color = floor(findgen(nfiles)/float(nfiles-1)*255)
+;for ifile=0,nfiles-1,nskip do begin
+tsave = [0.]
+toplot = [1,nfiles-1]
+for iplot=0,1 do begin
+ ifile = toplot[iplot]
+ inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
+ time = h5ra(inf, 'Header','Time')
+ tsave = [tsave, time]
+ print,' time= ',time
+ p = h5rd(inf,'PartType0/Coordinates')
+ v = h5rd(inf,'PartType0/Velocities')
+ id = h5rd(inf,'PartType0/ParticleIDs')
+ rho = h5rd(inf,'PartType0/Density')
+ h = h5rd(inf,'PartType0/SmoothingLength')
+ utherm = h5rd(inf,'PartType0/InternalEnergy')
+ indx = sort(id)
+
+; substract disk centre
+ for ic=0,2 do p[ic,*]=p[ic,*] - pcentre[ic]
+
+
+;; ; if you want to sort particles by ID
+;; id = id[indx]
+;; rho = rho[indx]
+;; utherm = utherm[indx]
+;; h = h[indx]
+;; for ic=0,2 do begin
+;; tmp = reform(p[ic,*]) & p[ic,*] = tmp[indx]
+;; tmp = reform(v[ic,*]) & v[ic,*] = tmp[indx]
+;; endfor
+
+ ip = floor(randomu(ifile+1,nplot)*n_elements(rho))
+ color = red
+ if(ifile eq 1) then begin
+ color=black
+ endif else begin
+ color=red
+ endelse
+ oplot,abs(p[2,ip]), rho[ip], psym=3, color=color
+
+ isave = isave + 1
+
+endfor
+
+; time in units of dynamical time
+tsave = tsave[1:*] / t_dyn
+
+label = ['']
+for i=0,n_elements(tsave)-1 do label=[label,'time/t_dynamic='+string(tsave[i],format='(f8.0)')]
+label = label[1:*]
+legend,['analytic',label[0],label[1]],linestyle=[0,0,0],color=[blue,black,red],box=0,/top,/right
+
+; make histograms of particle velocities
+xr = 1d-3 * [-1,1]
+bsize = 1.d-5
+ohist,v[0,*]/soundspeed,x,vx,xr[0],xr[1],bsize
+ohist,v[1,*]/soundspeed,y,vy,xr[0],xr[1],bsize
+ohist,v[2,*]/soundspeed,z,vz,xr[0],xr[1],bsize
+wset,2
+plot,x,vx,psym=10,xtitle='velocity/soundspeed',ytitle='pdf',/nodata,xr=xr,/xs
+oplot,x,vx,psym=10,color=black
+oplot,y,vy,psym=10,color=blue
+oplot,z,vz,psym=10,color=red
+legend,['vx/c','vy/c','vz/c'],linestyle=[0,0,0],color=[black,blue,red],box=0,/top,/right
+end
+
+
diff --git a/examples/Disk-Patch/HydroStatic/getGlass.sh b/examples/Disk-Patch/HydroStatic/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ffd92e88deae6e91237059adac2a6c2067caee46
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
diff --git a/examples/Disk-Patch/HydroStatic/makeIC.py b/examples/Disk-Patch/HydroStatic/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..7404f76aeaa485037b19052d84636c1656847a2b
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/makeIC.py
@@ -0,0 +1,254 @@
+###############################################################################
+ # 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
+import numpy
+import math
+import random
+import matplotlib.pyplot as plt
+
+# Generates a disk-patch in hydrostatic equilibrium
+# see Creasey, Theuns & Bower, 2013, for the equations:
+# disk parameters are: surface density sigma
+# scale height b
+# density: rho(z) = (sigma/2b) sech^2(z/b)
+# isothermal velocity dispersion = <v_z^2? = b pi G sigma
+# grad potential = 2 pi G sigma tanh(z/b)
+# potential = ln(cosh(z/b)) + const
+# Dynamical time = sqrt(b / (G sigma))
+# to obtain the 1/ch^2(z/b) profile from a uniform profile (a glass, say, or a uniform random variable), note that, when integrating in z
+# \int 0^z dz/ch^2(z) = tanh(z)-tanh(0) = \int_0^x dx = x (where the last integral refers to a uniform density distribution), so that z = atanh(x)
+# usage: python makeIC.py 1000
+
+# physical constants in cgs
+NEWTON_GRAVITY_CGS = 6.672e-8
+SOLAR_MASS_IN_CGS = 1.9885e33
+PARSEC_IN_CGS = 3.0856776e18
+PROTON_MASS_IN_CGS = 1.6726231e24
+YEAR_IN_CGS = 3.154e+7
+
+# 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)
+
+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
+
+
+# parameters of potential
+surface_density = 10.
+scale_height = 100.
+gamma = 5./3.
+
+# 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 'G=', const_G
+utherm = math.pi * const_G * surface_density * scale_height / (gamma-1)
+v_disp = numpy.sqrt(2 * utherm)
+soundspeed = numpy.sqrt(utherm / (gamma * (gamma-1.)))
+t_dyn = numpy.sqrt(scale_height / (const_G * surface_density))
+t_cross = scale_height / soundspeed
+print 'dynamical time = ',t_dyn,' sound crossing time = ',t_cross,' sound speed= ',soundspeed,' 3D velocity dispersion = ',v_disp,' thermal_energy= ',utherm
+
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 400. # [kpc]
+Radius = 100. # maximum radius of particles [kpc]
+G = const_G
+
+# File
+fileName = "Disk-Patch.hdf5"
+
+#---------------------------------------------------
+mass = 1
+
+#--------------------------------------------------
+
+
+# using glass ICs
+# read glass file and generate gas positions and tile it ntile times in each dimension
+ntile = 1
+inglass = 'glassCube_32.hdf5'
+infile = h5py.File(inglass, "r")
+one_glass_p = infile["/PartType0/Coordinates"][:,:]
+one_glass_h = infile["/PartType0/SmoothingLength"][:]
+
+# scale in [-0.5,0.5]*BoxSize / ntile
+one_glass_p[:,:] -= 0.5
+one_glass_p *= boxSize / ntile
+one_glass_h *= boxSize / ntile
+ndens_glass = (one_glass_h.shape[0]) / (boxSize/ntile)**3
+h_glass = numpy.amin(one_glass_h) * (boxSize/ntile)
+
+glass_p = []
+glass_h = []
+for ix in range(0,ntile):
+ for iy in range(0,ntile):
+ for iz in range(0,ntile):
+ shift = one_glass_p.copy()
+ shift[:,0] += (ix-(ntile-1)/2.) * boxSize / ntile
+ shift[:,1] += (iy-(ntile-1)/2.) * boxSize / ntile
+ shift[:,2] += (iz-(ntile-1)/2.) * boxSize / ntile
+ glass_p.append(shift)
+ glass_h.append(one_glass_h.copy())
+
+glass_p = numpy.concatenate(glass_p, axis=0)
+glass_h = numpy.concatenate(glass_h, axis=0)
+
+# random shuffle of glas ICs
+numpy.random.seed(12345)
+indx = numpy.random.rand(numpy.shape(glass_h)[0])
+indx = numpy.argsort(indx)
+glass_p = glass_p[indx, :]
+glass_h = glass_h[indx]
+
+# select numGas of them
+numGas = 8192
+pos = glass_p[0:numGas,:]
+h = glass_h[0:numGas]
+numGas = numpy.shape(pos)[0]
+
+# compute furthe properties of ICs
+column_density = surface_density * numpy.tanh(boxSize/2./scale_height)
+enclosed_mass = column_density * boxSize * boxSize
+pmass = enclosed_mass / numGas
+meanrho = enclosed_mass / boxSize**3
+print 'pmass= ',pmass,' mean(rho) = ', meanrho,' entropy= ', (gamma-1) * utherm / meanrho**(gamma-1)
+
+# desired density
+rho = surface_density / (2. * scale_height) / numpy.cosh(abs(pos[:,2])/scale_height)**2
+u = (1. + 0 * h) * utherm
+entropy = (gamma-1) * u / rho**(gamma-1)
+mass = 0.*h + pmass
+entropy_flag = 0
+vel = 0 + 0 * pos
+
+# move centre of disk to middle of box
+pos[:,:] += boxSize/2
+
+
+# create numPart dm particles
+numPart = 0
+
+# Create and write output file
+
+#File
+file = h5py.File(fileName, 'w')
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs
+grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] = [numGas, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numGas, numPart, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [entropy_flag]
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+
+# write gas particles
+grp0 = file.create_group("/PartType0")
+
+ds = grp0.create_dataset('Coordinates', (numGas, 3), 'f')
+ds[()] = pos
+
+ds = grp0.create_dataset('Velocities', (numGas, 3), 'f')
+ds[()] = vel
+
+ds = grp0.create_dataset('Masses', (numGas,), 'f')
+ds[()] = mass
+
+ds = grp0.create_dataset('SmoothingLength', (numGas,), 'f')
+ds[()] = h
+
+ds = grp0.create_dataset('InternalEnergy', (numGas,), 'f')
+u = numpy.full((numGas, ), utherm)
+if (entropy_flag == 1):
+ ds[()] = entropy
+else:
+ ds[()] = u
+
+ids = 1 + numpy.linspace(0, numGas, numGas, endpoint=False, dtype='L')
+ds = grp0.create_dataset('ParticleIDs', (numGas, ), 'L')
+ds[()] = ids
+
+print "Internal energy:", u[0]
+
+# generate dark matter particles if needed
+if(numPart > 0):
+
+ # set seed for random number
+ numpy.random.seed(1234)
+
+ grp1 = file.create_group("/PartType1")
+
+ radius = Radius * (numpy.random.rand(N))**(1./3.)
+ ctheta = -1. + 2 * numpy.random.rand(N)
+ stheta = numpy.sqrt(1.-ctheta**2)
+ phi = 2 * math.pi * numpy.random.rand(N)
+ r = numpy.zeros((numPart, 3))
+
+ speed = vrot
+ v = numpy.zeros((numPart, 3))
+ omega = speed / radius
+ period = 2.*math.pi/omega
+ print 'period = minimum = ',min(period), ' maximum = ',max(period)
+
+ v[:,0] = -omega * r[:,1]
+ v[:,1] = omega * r[:,0]
+
+ ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ ds[()] = r
+
+ ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+ ds[()] = v
+ v = numpy.zeros(1)
+
+ m = numpy.full((numPart, ),10)
+ ds = grp1.create_dataset('Masses', (numPart,), 'f')
+ ds[()] = m
+ m = numpy.zeros(1)
+
+ ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
+ ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
+ ds[()] = ids
+
+
+file.close()
+
+sys.exit()
diff --git a/examples/Disk-Patch/HydroStatic/test.pro b/examples/Disk-Patch/HydroStatic/test.pro
new file mode 100644
index 0000000000000000000000000000000000000000..31e027e3a308b04f1c59222e1a339786857061ac
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/test.pro
@@ -0,0 +1,142 @@
+;
+; 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 = 'Disk-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.,200.] * 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,'PartType0/ParticleIDs')
+nfollow = n_elements(id)
+
+; follow a subset
+; nfollow = min(4000, nfollow) ; 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
+vzout = fltarr(nfollow, nsave) ; z
+rout = fltarr(nfollow, nsave) ; rho
+hout = fltarr(nfollow, nsave) ; h
+uout = fltarr(nfollow, nsave) ; thermal energy
+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,'PartType0/Coordinates')
+ v = h5rd(inf,'PartType0/Velocities')
+ id = h5rd(inf,'PartType0/ParticleIDs')
+ rho = h5rd(inf,'PartType0/Density')
+ h = h5rd(inf,'PartType0/SmoothingLength')
+ utherm = h5rd(inf,'PartType0/InternalEnergy')
+ indx = sort(id)
+
+; if you want to sort particles by ID
+ id = id[indx]
+ rho = rho[indx]
+ utherm = utherm[indx]
+ h = h[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]
+ vzout[*,isave]= v[2,0:nfollow-1]
+ rout[*,isave] = rho[0:nfollow-1]
+ hout[*,isave] = h[0:nfollow-1]
+ uout[*,isave] = utherm[0:nfollow-1]
+ 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,*])
+
+
+; plot density profile and compare to analytic profile
+nplot = nfollow
+
+ ; plot density profile
+wset,0
+xr = [0, 3*scale_height]
+nbins = 100
+zpos = findgen(nbins)/float(nbins-1) * max(xr)
+dens = (surface_density/(2.d0*scale_height)) * 1./cosh(zpos/scale_height)^2
+plot,[0],[0],xr=xr,/xs,yr=[0,max(dens)*1.4],/ys,/nodata,xtitle='|z|',ytitle='density'
+oplot,zpos,dens,color=black,thick=3
+;oplot,abs(zout[*,1]),rout[*,1],psym=3 ; initial profile
+oplot,abs(zout[*,nsave-1]),rout[*,nsave-1],psym=3,color=red
+
+
+end
+
+
diff --git a/examples/IsothermalPotential/GravityOnly/test.pro b/examples/IsothermalPotential/GravityOnly/test.pro
index ac1f5e915c38c107408ffed9579c52944295f079..edfa50121d2e5adb7e039f3c38d6d4c0b4d5e34f 100644
--- a/examples/IsothermalPotential/GravityOnly/test.pro
+++ b/examples/IsothermalPotential/GravityOnly/test.pro
@@ -127,7 +127,7 @@ print,' relative Lz change: (per cent) ',minmax(dl) * 100.
if(iplot eq 1) then begin
; plot results on energy conservation for some particles
nplot = min(10, nfollow)
- wset,0
+ 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 (%)'
@@ -137,7 +137,7 @@ if(iplot eq 1) then begin
screen_to_png,'e-time.png'
; plot orbits of those particles
- wset,2
+ win,2
xr = [-100,100]
yr = xr
plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/iso,/nodata,xtitle='x',ytitle='y'
@@ -146,7 +146,7 @@ if(iplot eq 1) then begin
screen_to_png,'orbit.png'
; plot radial position of these particles
- wset,4
+ win,4
xr = [min(tout), max(tout)]
yr = [0,80]
plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='t',ytitle='r'
@@ -155,7 +155,7 @@ for i=0,nplot-1 do begin dr = sqrt(reform(xout[i,*])^2 + reform(yout[i,*])^2) &
screen_to_png,'r-time.png'
; make histogram of energy changes at end
- wset,6
+ 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'
diff --git a/examples/main.c b/examples/main.c
index efd40c3da9802ebfab72a51eb841c035c8490eba..c9253ca4dd40ce787b62b0c8ebf89eebbe08cb2c 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -332,7 +332,8 @@ int main(int argc, char *argv[]) {
/* Initialise the external potential properties */
struct external_potential potential;
- if (with_external_gravity) potential_init(params, &us, &potential);
+ if (with_external_gravity)
+ potential_init(params, &prog_const, &us, &potential);
if (with_external_gravity && myrank == 0) potential_print(&potential);
/* Read particles and space information from (GADGET) ICs */
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index edb0885d621975850a09e4298b8e035ebb45a3cd..0a15f5e1628f70306cf021c4bd05bef699d78214 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -66,15 +66,22 @@ DomainDecomposition:
# Point mass external potentials
PointMass:
- position_x: 50. # location of external point mass in internal units
+ position_x: 50. # location of external point mass (internal units)
position_y: 50.
position_z: 50.
- mass: 1e10 # mass of external point mass in internal units
- timestep_mult: 0.03 # Pre-factor for the time-step condition
+ mass: 1e10 # mass of external point mass (internal units)
+ timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
IsothermalPotential:
- position_x: 100. # Location of centre of isothermal potential in internal units
+ position_x: 100. # Location of centre of isothermal potential (internal units)
position_y: 100.
position_z: 100.
- vrot: 200. # Rotation speed of isothermal potential in internal units
- timestep_mult: 0.03 # Pre-factor for the time-step condition
+ vrot: 200. # Rotation speed of isothermal potential (internal units)
+ timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
+
+Disk-PatchPotential:
+ surface_density: 10. # Surface density of the disk (internal units)
+ scale_height: 100. # Scale height of the disk (internal units)
+ z_disk: 200. # Disk height (internal units)
+ timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
+ growth_time: 5. # (Optional) Time for the disk to grow to its final size (multiple of the dynamical time)
diff --git a/src/const.h b/src/const.h
index 6710fe9e52dfdbede468282d84d4441b923a59fe..cbaa0054fe3bc9f829c59b19b71197264fd0c5af 100644
--- a/src/const.h
+++ b/src/const.h
@@ -36,6 +36,9 @@
/* Time integration constants. */
#define const_max_u_change 0.1f
+/* Thermal energy per unit mass used as a constant for the isothermal EoS */
+#define const_isothermal_internal_energy 20.2615290634f
+
/* Dimensionality of the problem */
#define HYDRO_DIMENSION_3D
//#define HYDRO_DIMENSION_2D
@@ -72,6 +75,7 @@
/* External gravity properties */
#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+//#define EXTERNAL_POTENTIAL_DISK_PATCH
/* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS
diff --git a/src/equation_of_state.h b/src/equation_of_state.h
index b4a36e8a3ef0bcc281d1f939f89fde08ecf00be9..af59d8a2cad1632c67b6d377b5ed9dfe9484b4aa 100644
--- a/src/equation_of_state.h
+++ b/src/equation_of_state.h
@@ -132,66 +132,93 @@ gas_soundspeed_from_internal_energy(float density, float u) {
/**
* @brief Returns the internal energy given density and entropy
*
- * @param density The density
- * @param entropy The entropy
+ * Since we are using an isothermal EoS, the entropy value is ignored
+ * Computes \f$u = u_{cst}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
*/
__attribute__((always_inline)) INLINE static float
gas_internal_energy_from_entropy(float density, float entropy) {
- error("Missing definition !");
- return 0.f;
+ return const_isothermal_internal_energy;
}
-
/**
* @brief Returns the pressure given density and entropy
*
- * @param density The density
- * @param entropy The entropy
+ * Since we are using an isothermal EoS, the entropy value is ignored
+ * Computes \f$P = (\gamma - 1)u_{cst}\rho\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
*/
__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
float density, float entropy) {
- error("Missing definition !");
- return 0.f;
+ return hydro_gamma_minus_one * const_isothermal_internal_energy * density;
+}
+
+/**
+ * @brief Returns the sound speed given density and entropy
+ *
+ * Since we are using an isothermal EoS, the entropy value is ignored
+ * Computes \f$c = \sqrt{u_{cst} \gamma \rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
+ float density, float entropy) {
+
+ return sqrtf(const_isothermal_internal_energy * hydro_gamma *
+ hydro_gamma_minus_one);
}
/**
* @brief Returns the entropy given density and internal energy
*
- * @param density The density
- * @param u The internal energy
+ * Since we are using an isothermal EoS, the energy value is ignored
+ * Computes \f$S = \frac{(\gamma - 1)u_{cst}}{\rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
*/
__attribute__((always_inline)) INLINE static float
gas_entropy_from_internal_energy(float density, float u) {
- error("Missing definition !");
- return 0.f;
+ return hydro_gamma_minus_one * const_isothermal_internal_energy *
+ pow_minus_gamma_minus_one(density);
}
/**
* @brief Returns the pressure given density and internal energy
*
- * @param density The density
- * @param u The internal energy
+ * Since we are using an isothermal EoS, the energy value is ignored
+ * Computes \f$P = (\gamma - 1)u_{cst}\rho\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
*/
__attribute__((always_inline)) INLINE static float
gas_pressure_from_internal_energy(float density, float u) {
- error("Missing definition !");
- return 0.f;
+ return hydro_gamma_minus_one * const_isothermal_internal_energy * density;
}
/**
* @brief Returns the sound speed given density and internal energy
*
- * @param density The density
- * @param u The internal energy
+ * Since we are using an isothermal EoS, the energy value is ignored
+ * Computes \f$c = \sqrt{u_{cst} \gamma \rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
*/
__attribute__((always_inline)) INLINE static float
gas_soundspeed_from_internal_energy(float density, float u) {
- error("Missing definition !");
- return 0.f;
+ return sqrtf(const_isothermal_internal_energy * hydro_gamma *
+ hydro_gamma_minus_one);
}
/* ------------------------------------------------------------------------- */
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index f9b67c96331e7572b0200093f0c32eee5d2391cd..da6b9aa01d7d382b1d6aa85f732f1fae1169ac34 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -48,7 +48,10 @@ gravity_compute_timestep_external(const struct external_potential* potential,
dt = fminf(dt, external_gravity_isothermalpotential_timestep(potential,
phys_const, gp));
#endif
-
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+ dt = fminf(dt,
+ external_gravity_disk_patch_timestep(potential, phys_const, gp));
+#endif
return dt;
}
@@ -128,12 +131,13 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
*
* This function only branches towards the potential chosen by the user.
*
+ * @param time The current time in internal units.
* @param potential The properties of the external potential.
* @param phys_const The physical constants in internal units.
* @param gp The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void external_gravity(
- const struct external_potential* potential,
+ double time, const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* gp) {
#ifdef EXTERNAL_POTENTIAL_POINTMASS
@@ -142,6 +146,9 @@ __attribute__((always_inline)) INLINE static void external_gravity(
#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
external_gravity_isothermalpotential(potential, phys_const, gp);
#endif
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+ external_gravity_disk_patch_potential(time, potential, phys_const, gp);
+#endif
}
/**
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 0e4eaf150c764e27156747bb01db20a71f03c7b5..815969975b4f5e6b39099e71bbbec4e43c875ddc 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -56,6 +56,15 @@ void hydro_props_init(struct hydro_props *p,
void hydro_props_print(const struct hydro_props *p) {
+#if defined(EOS_IDEAL_GAS)
+ message("Equation of state: Ideal gas.");
+#elif defined(EOS_ISOTHERMAL_GAS)
+ message(
+ "Equation of state: Isothermal with internal energy "
+ "per unit mass set to %f.",
+ const_isothermal_internal_energy);
+#endif
+
message("Adiabatic index gamma: %f.", hydro_gamma);
message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
diff --git a/src/potentials.c b/src/potentials.c
index 5cbe05e008b7de17310b1e738e032af90684c25e..dd7aed8712c01921a01462bf9de713433c81f5c1 100644
--- a/src/potentials.c
+++ b/src/potentials.c
@@ -29,11 +29,13 @@
* of units.
*
* @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
* @param us The current internal system of units
* @param potential The external potential properties to initialize
*/
void potential_init(const struct swift_params* parameter_file,
- struct UnitSystem* us,
+ const struct phys_const* phys_const,
+ const struct UnitSystem* us,
struct external_potential* potential) {
#ifdef EXTERNAL_POTENTIAL_POINTMASS
@@ -65,6 +67,22 @@ void potential_init(const struct swift_params* parameter_file,
parameter_file, "IsothermalPotential:timestep_mult");
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+ potential->disk_patch_potential.surface_density = parser_get_param_double(
+ parameter_file, "Disk-PatchPotential:surface_density");
+ potential->disk_patch_potential.scale_height = parser_get_param_double(
+ parameter_file, "Disk-PatchPotential:scale_height");
+ potential->disk_patch_potential.z_disk =
+ parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk");
+ potential->disk_patch_potential.timestep_mult = parser_get_param_double(
+ parameter_file, "Disk-PatchPotential:timestep_mult");
+ potential->disk_patch_potential.growth_time = parser_get_opt_param_double(
+ parameter_file, "Disk-PatchPotential:growth_time", 0.);
+ potential->disk_patch_potential.dynamical_time =
+ sqrt(potential->disk_patch_potential.scale_height /
+ (phys_const->const_newton_G *
+ potential->disk_patch_potential.surface_density));
+#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
}
/**
@@ -78,7 +96,7 @@ void potential_print(const struct external_potential* potential) {
message(
"Point mass properties are (x,y,z) = (%e, %e, %e), M = %e timestep "
- "multiplier = %e",
+ "multiplier = %e.",
potential->point_mass.x, potential->point_mass.y, potential->point_mass.z,
potential->point_mass.mass, potential->point_mass.timestep_mult);
@@ -88,10 +106,25 @@ void potential_print(const struct external_potential* potential) {
message(
"Isothermal potential properties are (x,y,z) = (%e, %e, %e), vrot = %e "
- "timestep multiplier= %e",
+ "timestep multiplier = %e.",
potential->isothermal_potential.x, potential->isothermal_potential.y,
potential->isothermal_potential.z, potential->isothermal_potential.vrot,
potential->isothermal_potential.timestep_mult);
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
+
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+
+ message(
+ "Disk-patch potential properties are surface_density = %e disk height= "
+ "%e scale height = %e timestep multiplier = %e.",
+ potential->disk_patch_potential.surface_density,
+ potential->disk_patch_potential.z_disk,
+ potential->disk_patch_potential.scale_height,
+ potential->disk_patch_potential.timestep_mult);
+
+ if (potential->disk_patch_potential.growth_time > 0.)
+ message("Disk will grow for %f dynamiccal times.",
+ potential->disk_patch_potential.growth_time);
+#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
}
diff --git a/src/potentials.h b/src/potentials.h
index 5373cc1f3f55b0c6d9876cbe6348fdcefbe242aa..f0de5b4dc784be2e613c634bf687300d57f40ecd 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -25,6 +25,7 @@
#include "../config.h"
/* Some standard headers. */
+#include <float.h>
#include <math.h>
/* Local includes. */
@@ -53,13 +54,123 @@ struct external_potential {
float timestep_mult;
} isothermal_potential;
#endif
+
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+ struct {
+ double surface_density;
+ double scale_height;
+ double z_disk;
+ double dynamical_time;
+ double growth_time;
+ double timestep_mult;
+ } disk_patch_potential;
+#endif
};
-/* Include external isothermal potential */
+/* ------------------------------------------------------------------------- */
+
+/* external potential: disk_patch */
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+
+/**
+ * @brief Computes the time-step from the acceleration due to a hydrostatic
+ * disk.
+ *
+ * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
+ *
+ * @param phys_cont The physical constants in internal units.
+ * @param gp Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+external_gravity_disk_patch_timestep(const struct external_potential* potential,
+ const struct phys_const* const phys_const,
+ const struct gpart* const g) {
+
+ /* initilize time step to disk dynamical time */
+ const float dt_dyn = potential->disk_patch_potential.dynamical_time;
+ float dt = dt_dyn;
+
+ /* absolute value of height above disk */
+ const float dz = fabs(g->x[2] - potential->disk_patch_potential.z_disk);
+
+ /* vertical cceleration */
+ const float z_accel = 2 * M_PI * phys_const->const_newton_G *
+ potential->disk_patch_potential.surface_density *
+ tanh(dz / potential->disk_patch_potential.scale_height);
+
+ /* demand that dt * velocity < fraction of scale height of disk */
+ float dt1 = FLT_MAX;
+ if (fabs(g->v_full[2]) > 0) {
+ dt1 = potential->disk_patch_potential.scale_height / fabs(g->v_full[2]);
+ if (dt1 < dt) dt = dt1;
+ }
+
+ /* demand that dt^2 * acceleration < fraction of scale height of disk */
+ float dt2 = FLT_MAX;
+ if (fabs(z_accel) > 0) {
+ dt2 = potential->disk_patch_potential.scale_height / fabs(z_accel);
+ if (dt2 < dt * dt) dt = sqrt(dt2);
+ }
+
+ /* demand that dt^3 jerk < fraction of scale height of disk */
+ float dt3 = FLT_MAX;
+ if (abs(g->v_full[2]) > 0) {
+ const float dz_accel_over_dt =
+ 2 * M_PI * phys_const->const_newton_G *
+ potential->disk_patch_potential.surface_density /
+ potential->disk_patch_potential.scale_height /
+ cosh(dz / potential->disk_patch_potential.scale_height) /
+ cosh(dz / potential->disk_patch_potential.scale_height) *
+ fabs(g->v_full[2]);
+
+ dt3 = potential->disk_patch_potential.scale_height / fabs(dz_accel_over_dt);
+ if (dt3 < dt * dt * dt) dt = pow(dt3, 1. / 3.);
+ }
+
+ return potential->disk_patch_potential.timestep_mult * dt;
+}
+
+/**
+ * @brief Computes the gravitational acceleration along z due to a hydrostatic
+ * disk
+ *
+ * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
+ *
+ * @param phys_cont The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static void
+external_gravity_disk_patch_potential(
+ const double time, const struct external_potential* potential,
+ const struct phys_const* const phys_const, struct gpart* g) {
+
+ const float G_newton = phys_const->const_newton_G;
+ const float dz = g->x[2] - potential->disk_patch_potential.z_disk;
+ const float t_dyn = potential->disk_patch_potential.dynamical_time;
+
+ float reduction_factor = 1.;
+ if (time < potential->disk_patch_potential.growth_time * t_dyn)
+ reduction_factor =
+ time / (potential->disk_patch_potential.growth_time * t_dyn);
+
+ const float z_accel =
+ reduction_factor * 2 * G_newton * M_PI *
+ potential->disk_patch_potential.surface_density *
+ tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
+
+ /* Accelerations. Note that they are multiplied by G later on */
+ if (dz > 0) g->a_grav[2] -= z_accel / G_newton;
+ if (dz < 0) g->a_grav[2] += z_accel / G_newton;
+}
+#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
+
+/* ------------------------------------------------------------------------- */
+
#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
/**
- * @brief Computes the time-step due to the acceleration from a point mass
+ * @brief Computes the time-step due to the acceleration from an isothermal
+ * potential.
*
* @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units.
@@ -73,6 +184,7 @@ external_gravity_isothermalpotential_timestep(
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
+
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const float drdv =
dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
@@ -92,10 +204,10 @@ external_gravity_isothermalpotential_timestep(
}
/**
- * @brief Computes the gravitational acceleration of a particle due to a point
- * mass
+ * @brief Computes the gravitational acceleration from an isothermal potential.
*
- * Note that the accelerations are multiplied by Newton's G constant later on.
+ * Note that the accelerations are multiplied by Newton's G constant
+ * later on.
*
* @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units.
@@ -107,7 +219,6 @@ external_gravity_isothermalpotential(const struct external_potential* potential,
struct gpart* g) {
const float G_newton = phys_const->const_newton_G;
-
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
@@ -119,11 +230,12 @@ external_gravity_isothermalpotential(const struct external_potential* potential,
g->a_grav[0] += term * dx;
g->a_grav[1] += term * dy;
g->a_grav[2] += term * dz;
- // error(" %f %f %f %f", vrot, rinv2, dx, g->a_grav[0]);
}
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
+/* ------------------------------------------------------------------------- */
+
/* Include exteral pointmass potential */
#ifdef EXTERNAL_POTENTIAL_POINTMASS
@@ -161,10 +273,11 @@ external_gravity_pointmass_timestep(const struct external_potential* potential,
}
/**
- * @brief Computes the gravitational acceleration of a particle due to a point
- * mass
+ * @brief Computes the gravitational acceleration of a particle due to a
+ * point mass
*
- * Note that the accelerations are multiplied by Newton's G constant later on.
+ * Note that the accelerations are multiplied by Newton's G constant later
+ * on.
*
* @param potential The proerties of the external potential.
* @param phys_const The physical constants in internal units.
@@ -187,9 +300,12 @@ __attribute__((always_inline)) INLINE static void external_gravity_pointmass(
#endif /* EXTERNAL_POTENTIAL_POINTMASS */
+/* ------------------------------------------------------------------------- */
+
/* Now, some generic functions, defined in the source file */
void potential_init(const struct swift_params* parameter_file,
- struct UnitSystem* us,
+ const struct phys_const* phys_const,
+ const struct UnitSystem* us,
struct external_potential* potential);
void potential_print(const struct external_potential* potential);
diff --git a/src/runner.c b/src/runner.c
index 23fd820d119230a499a87a447562fc631142655e..eabc622e8c3aa41b5b06c604a4a540b8cf7c13b9 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -105,7 +105,7 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
const int ti_current = r->e->ti_current;
const struct external_potential *potential = r->e->external_potential;
const struct phys_const *constants = r->e->physical_constants;
-
+ const double time = r->e->time;
TIMER_TIC;
/* Recurse? */
@@ -128,7 +128,7 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
/* Is this part within the time step? */
if (g->ti_end <= ti_current) {
- external_gravity(potential, constants, g);
+ external_gravity(time, potential, constants, g);
}
}