Skip to content
Snippets Groups Projects
Commit 2d5efd2a authored by Tom Theuns's avatar Tom Theuns
Browse files

updated dis-patch. Verified stable equilibirum over 250 dynamical times with 8192 particles

parent 668bb3ce
Branches
Tags
1 merge request!217Disk patch
# 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
h_scaling: 1. # A scaling factor to apply to all smoothing lengths in the ICs.
shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
shift_z: 0.
# External potential parameters
PointMass:
position_x: 50. # location of external point mass in internal units
position_y: 50.
position_z: 50.
mass: 1e10 # mass of external point mass in internal units
timestep_mult: 0.03 # controls time step
IsothermalPotential:
position_x: 100. # location of centre of isothermal potential in internal units
position_y: 100.
position_z: 100.
vrot: 200. # rotation speed of isothermal potential in internal units
timestep_mult: 0.03 # controls time step
Disk-PatchPotential:
surface_density: 10.
scale_height: 100.
z_disk: 200.
timestep_mult: 0.03
...@@ -19,9 +19,14 @@ uV = 1d5 ; unit of velocity ...@@ -19,9 +19,14 @@ uV = 1d5 ; unit of velocity
surface_density = 10. surface_density = 10.
scale_height = 100. scale_height = 100.
z_disk = 200.; z_disk = 200.;
gamma = 5./3.
; derived units ; derived units
constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ; constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
pcentre = [0.,0.,z_disk] * pc / 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 + '*' infile = indir + basefile + '*'
...@@ -48,7 +53,7 @@ rbins = (surface_density/(2.*scale_height)) / cosh(abs(zbins)/scale_height)^2 ...@@ -48,7 +53,7 @@ rbins = (surface_density/(2.*scale_height)) / cosh(abs(zbins)/scale_height)^2
; plot analytic profile ; plot analytic profile
wset,0 wset,0
plot,[0],[0],xr=[0,2*scale_height],yr=[0,max(rbins)],/nodata,xtitle='|z|',ytitle=textoidl('\rho') plot,[0],[0],xr=[0,2*scale_height],yr=[0,max(rbins)],/nodata,xtitle='|z|',ytitle=textoidl('\rho')
oplot,zbins,rbins oplot,zbins,rbins,color=blue
ifile = 0 ifile = 0
nskip = nfiles - 1 nskip = nfiles - 1
...@@ -56,8 +61,10 @@ isave = 0 ...@@ -56,8 +61,10 @@ isave = 0
nplot = 8192 ; randomly plot particles nplot = 8192 ; randomly plot particles
color = floor(findgen(nfiles)/float(nfiles-1)*255) color = floor(findgen(nfiles)/float(nfiles-1)*255)
;for ifile=0,nfiles-1,nskip do begin ;for ifile=0,nfiles-1,nskip do begin
tsave = [0.] tsave = [0.]
for ifile=0,nfiles-1,nskip do begin toplot = [1,nfiles-1]
for iplot=0,1 do begin
ifile = toplot[iplot]
inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5' inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
time = h5ra(inf, 'Header','Time') time = h5ra(inf, 'Header','Time')
tsave = [tsave, time] tsave = [tsave, time]
...@@ -86,19 +93,37 @@ for ifile=0,nfiles-1,nskip do begin ...@@ -86,19 +93,37 @@ for ifile=0,nfiles-1,nskip do begin
ip = floor(randomu(ifile+1,nplot)*n_elements(rho)) ip = floor(randomu(ifile+1,nplot)*n_elements(rho))
color = red color = red
if(ifile eq 0) then color=black 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 oplot,abs(p[2,ip]), rho[ip], psym=3, color=color
isave = isave + 1 isave = isave + 1
endfor endfor
tsave = tsave[1:*]
; time in units of dynamical time
tsave = tsave[1:*] / t_dyn
label = [''] label = ['']
for i=0,n_elements(tsave)-1 do label=[label,'t='+string(tsave[i],format='(f8.0)')] for i=0,n_elements(tsave)-1 do label=[label,'time/t_dynamic='+string(tsave[i],format='(f8.0)')]
label = label[1:*] label = label[1:*]
legend,[label[0],label[1]],linestyle=[0,0],color=[black,red],box=0,/top,/right 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 end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment