Commit 602abd1b authored by Matthieu Schaller's avatar Matthieu Schaller

Clean up the gravity-only disc-patch example

parent 419a58fa
# 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
......@@ -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
......
......@@ -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
;
; 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
###############################################################################
# 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, "%")
......@@ -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)
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment