Commit 33787d31 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Removed the IDL scripts for the disc patch test case. Updated the python...

Removed the IDL scripts for the disc patch test case. Updated the python script to use the new range of snapshots.
parent 8872a1f5
......@@ -18,3 +18,5 @@ output to 'Disc-Patch-dynamic.hdf5'. These are now the ICs for the actual test.
When running SWIFT with the parameters from 'disc-patch.yml' and an
ideal gas EoS on these ICs the disc should stay in equilibrium.
The solution can be checked using the 'plotSolution.py' script.
;
; 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-dynamic_'
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 = 100. ; surface density of all mass, which generates the gravitational potential
scale_height = 100.
z_disk = 200. ;
fgas = 0.1 ; gas fraction
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.)
temp = (utherm*uV^2)*phys.m_h/phys.kb
soundspeed = sqrt(gamma * (gamma-1.) * utherm)
t_dyn = sqrt(scale_height / (constG * surface_density))
rho0 = fgas*(surface_density)/(2.*scale_height)
print,' dynamical time = ',t_dyn,' = ',t_dyn*UL/uV/(1d6*phys.yr),' Myr'
print,' thermal energy per unit mass = ',utherm
print,' central density = ',rho0,' = ',rho0*uM/uL^3/m_h,' particles/cm^3'
print,' central temperature = ',temp
lambda = 2 * !pi * phys.G^1.5 * (scale_height*uL)^1.5 * (surface_density * uM/uL^2)^0.5 * phys.m_h^2 / (gamma-1) / fgas
print,' lambda = ',lambda
stop
;
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
......@@ -20,7 +20,7 @@
##
# This script plots the Disc-Patch_*.hdf5 snapshots.
# It takes two (optional) parameters: the counter value of the first and last
# snapshot to plot (default: 0 81).
# snapshot to plot (default: 0 21).
##
import numpy as np
......@@ -34,12 +34,12 @@ import sys
# Parameters
surface_density = 10.
scale_height = 100.
z_disc = 200.
utherm = 20.2615290634
z_disc = 400.
utherm = 20.2678457288
gamma = 5. / 3.
start = 0
stop = 81
stop = 21
if len(sys.argv) > 1:
start = int(sys.argv[1])
if len(sys.argv) > 2:
......
;
; 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.,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
Supports Markdown
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