From 602abd1bb09e9fa4efbe868674b3d50cc0414c52 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Sun, 12 Apr 2020 17:08:41 +0200
Subject: [PATCH] Clean up the gravity-only disc-patch example

---
 .../DiscPatch/GravityOnly/disc-patch.yml      |  22 +--
 .../DiscPatch/GravityOnly/makeIC.py           |  70 ++++----
 .../GravityTests/DiscPatch/GravityOnly/run.sh |   4 +
 .../DiscPatch/GravityOnly/test.pro            | 158 ------------------
 .../DiscPatch/GravityOnly/test.py             |  92 ++++++++++
 examples/parameter_example.yml                |   8 +-
 6 files changed, 139 insertions(+), 215 deletions(-)
 delete mode 100644 examples/GravityTests/DiscPatch/GravityOnly/test.pro
 create mode 100644 examples/GravityTests/DiscPatch/GravityOnly/test.py

diff --git a/examples/GravityTests/DiscPatch/GravityOnly/disc-patch.yml b/examples/GravityTests/DiscPatch/GravityOnly/disc-patch.yml
index bcc7d1a3d..34dcfdf4a 100644
--- a/examples/GravityTests/DiscPatch/GravityOnly/disc-patch.yml
+++ b/examples/GravityTests/DiscPatch/GravityOnly/disc-patch.yml
@@ -1,8 +1,8 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.98848e33    # M_sun in grams
-  UnitLength_in_cgs:   3.08567758e18 # parsec in centimeters
-  UnitVelocity_in_cgs: 1e5           # km/s in cm/s
+  UnitMass_in_cgs:     1.98841e33       # M_sun in grams
+  UnitLength_in_cgs:   3.08567758149e18 # parsec in centimeters
+  UnitVelocity_in_cgs: 1e5              # km/s in cm/s
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 
@@ -23,14 +23,6 @@ Snapshots:
   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:  Disc-Patch.hdf5       # The file to read
@@ -40,5 +32,9 @@ InitialConditions:
 DiscPatchPotential:
   surface_density: 10.
   scale_height:    100.
-  z_disc:          300.
-  timestep_mult:   0.03
+  x_disc:          300.
+  timestep_mult:   0.01
+
+Scheduler:
+  max_top_level_cells:  8
+  tasks_per_cell:       5
diff --git a/examples/GravityTests/DiscPatch/GravityOnly/makeIC.py b/examples/GravityTests/DiscPatch/GravityOnly/makeIC.py
index 3abf4f87f..baeb2f359 100644
--- a/examples/GravityTests/DiscPatch/GravityOnly/makeIC.py
+++ b/examples/GravityTests/DiscPatch/GravityOnly/makeIC.py
@@ -19,29 +19,32 @@
  ##############################################################################
 
 import h5py
-import sys
 import numpy
+import sys
 import math
-import random
 
-# Generates N particles in a box of [0:BoxSize,0:BoxSize,-2scale_height:2scale_height]
+# Generates N particles in a box of [-2scale_height:2scale_height, 0:BoxSize, 0:BoxSize]
+#
+# The potential is long the x-axis. 
+#
 # see Creasey, Theuns & Bower, 2013, for the equations:
 # disc 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
+# density: rho(x) = (sigma/2b) sech^2(x/b)
+# isothermal velocity dispersion = <v_x^2? = b pi G sigma
 # grad potential  = 2 pi G sigma tanh(z/b)
-# potential       = ln(cosh(z/b)) + const
+# potential       = ln(cosh(x/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)
+# to obtain the 1/ch^2(x/b) profile from a uniform profile (a glass, say, or a uniform random variable), note that, when integrating in x
+# \int 0^x dx/ch^2(x) = tanh(x)-tanh(0) = \int_0^u du = u (where the last integral refers to a uniform density distribution), so that x = atanh(u)
+#
 # usage: python makeIC.py 1000 
 
 # physical constants in cgs
-NEWTON_GRAVITY_CGS  = 6.67408e-8
-SOLAR_MASS_IN_CGS   = 1.98848e33
-PARSEC_IN_CGS       = 3.08567758e18
-YEAR_IN_CGS         = 3.15569252e7
+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)
@@ -67,25 +70,14 @@ 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
+periodic = 1             # 1 For periodic box
+boxSize  = 600.          #  
+Radius   = 100.          # maximum radius of particles [kpc]
+G        = const_G 
+mass     = 1
+numPart  = int(sys.argv[1])  # Number of particles
 fileName = "Disc-Patch.hdf5" 
 
-
-#---------------------------------------------------
-numPart        = N
-mass           = 1
-internalEnergy = P / ((gamma - 1.)*rho)
-
 #--------------------------------------------------
 
 #File
@@ -115,25 +107,23 @@ grp.attrs["Dimension"] = 3
 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
+x      = scale_height * numpy.arctanh(numpy.random.rand(2*numPart))
+gd     = x < boxSize / 2
+r[:,0] = x[gd][0:numPart]
+random = numpy.random.rand(numPart) > 0.5
+r[random,0] *= -1
+r[:,0] += 0.5 * boxSize
+
+r[:,1] = numpy.random.rand(numPart) * boxSize
+r[:,2] = numpy.random.rand(numPart) * 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
diff --git a/examples/GravityTests/DiscPatch/GravityOnly/run.sh b/examples/GravityTests/DiscPatch/GravityOnly/run.sh
index ab05b6032..3b7c51683 100755
--- a/examples/GravityTests/DiscPatch/GravityOnly/run.sh
+++ b/examples/GravityTests/DiscPatch/GravityOnly/run.sh
@@ -7,4 +7,8 @@ then
     python makeIC.py 1000
 fi
 
+# Run SWIFT
 ../../../swift --external-gravity --threads=2 disc-patch.yml
+
+# Verify energy conservation
+python3 test.py
diff --git a/examples/GravityTests/DiscPatch/GravityOnly/test.pro b/examples/GravityTests/DiscPatch/GravityOnly/test.pro
deleted file mode 100644
index 04e0afdf7..000000000
--- a/examples/GravityTests/DiscPatch/GravityOnly/test.pro
+++ /dev/null
@@ -1,158 +0,0 @@
-;
-;  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
-
-
diff --git a/examples/GravityTests/DiscPatch/GravityOnly/test.py b/examples/GravityTests/DiscPatch/GravityOnly/test.py
new file mode 100644
index 000000000..e15047338
--- /dev/null
+++ b/examples/GravityTests/DiscPatch/GravityOnly/test.py
@@ -0,0 +1,92 @@
+###############################################################################
+# 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, "%")
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 831619622..962c806b6 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -280,13 +280,13 @@ NFWPotential:
   critical_density:   127.4    # Critical density (internal units).
   timestep_mult:      0.01     # Dimensionless pre-factor for the time-step condition, basically determines fraction of orbital time we need to do an integration step
 
-# Disk-patch potential parameters
+# Disk-patch potential parameters. The potential is along the x-axis.
 DiscPatchPotential:
   surface_density: 10.      # Surface density of the disc (internal units)
   scale_height:    100.     # Scale height of the disc (internal units)
-  z_disc:          400.     # Position of the disc along the z-axis (internal units)
-  z_trunc:         300.     # (Optional) Distance from the disc along z-axis above which the potential gets truncated.
-  z_max:           380.     # (Optional) Distance from the disc along z-axis above which the potential is set to 0.
+  x_disc:          400.     # Position of the disc along the z-axis (internal units)
+  x_trunc:         300.     # (Optional) Distance from the disc along z-axis above which the potential gets truncated.
+  x_max:           380.     # (Optional) Distance from the disc along z-axis above which the potential is set to 0.
   timestep_mult:   0.03     # Dimensionless pre-factor for the time-step condition
   growth_time:     5.       # (Optional) Time for the disc to grow to its final size (multiple of the dynamical time)
 
-- 
2.26.2