Skip to content
Snippets Groups Projects
Commit ad8aa1fb authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'disc_patch_revived' into 'master' 

Disc patch revived

Implements the changes discussed in #311. 

 - The disc patch example is now 400x400x800 in length with correct internal energy.
 - The accelerations get truncated at |z-z_disc| > 300 with a cosine transition
 - The accelerations are 0 at |z-z_disc| > 380.
 - The calculation is significantly faster thanks to the extraction of constant terms.

See merge request !349
parents 356a3348 1334f2e2
Branches
Tags
1 merge request!349Disc patch revived
Hide whitespace changes
Inline Side-by-side
Showing
with 362 additions and 548 deletions
examples/DiscPatch/HydroStatic/README
+ 2
− 0
View file @ ad8aa1fb
...@@ -18,3 +18,5 @@ output to 'Disc-Patch-dynamic.hdf5'. These are now the ICs for the actual test. ...@@ -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 When running SWIFT with the parameters from 'disc-patch.yml' and an
ideal gas EoS on these ICs the disc should stay in equilibrium. ideal gas EoS on these ICs the disc should stay in equilibrium.
The solution can be checked using the 'plotSolution.py' script.
examples/DiscPatch/HydroStatic/disc-patch-icc.yml
+ 12
− 10
View file @ ad8aa1fb
# Define the system of units to use internally. # Define the system of units to use internally.
InternalUnitSystem: InternalUnitSystem:
UnitMass_in_cgs: 1.9885e33 # Grams UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.0856776e18 # Centimeters UnitLength_in_cgs: 3.08567758149e18 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin UnitTemp_in_cgs: 1 # Kelvin
...@@ -11,17 +11,17 @@ TimeIntegration: ...@@ -11,17 +11,17 @@ TimeIntegration:
time_begin: 0 # The starting time of the simulation (in internal units). time_begin: 0 # The starting time of the simulation (in internal units).
time_end: 968. # The end 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_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). dt_max: 10. # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics # Parameters governing the conserved quantities statistics
Statistics: Statistics:
delta_time: 1 # Time between statistics output delta_time: 12. # Time between statistics output
# Parameters governing the snapshots # Parameters governing the snapshots
Snapshots: Snapshots:
basename: Disc-Patch # Common part of the name of output files basename: Disc-Patch # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units) time_first: 0. # Time of the first output (in internal units)
delta_time: 12. # Time difference between consecutive outputs (in internal units) delta_time: 48. # Time difference between outputs (in internal units)
# Parameters for the hydrodynamics scheme # Parameters for the hydrodynamics scheme
SPH: SPH:
...@@ -29,7 +29,7 @@ SPH: ...@@ -29,7 +29,7 @@ SPH:
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours. delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. 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_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). h_max: 60. # Maximal smoothing length allowed (in internal units).
# Parameters related to the initial conditions # Parameters related to the initial conditions
InitialConditions: InitialConditions:
...@@ -39,6 +39,8 @@ InitialConditions: ...@@ -39,6 +39,8 @@ InitialConditions:
DiscPatchPotential: DiscPatchPotential:
surface_density: 10. surface_density: 10.
scale_height: 100. scale_height: 100.
z_disc: 200. z_disc: 400.
z_trunc: 300.
z_max: 350.
timestep_mult: 0.03 timestep_mult: 0.03
growth_time: 5. growth_time: 5.
examples/DiscPatch/HydroStatic/disc-patch.yml
+ 12
− 10
View file @ ad8aa1fb
# Define the system of units to use internally. # Define the system of units to use internally.
InternalUnitSystem: InternalUnitSystem:
UnitMass_in_cgs: 1.9885e33 # Grams UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.0856776e18 # Centimeters UnitLength_in_cgs: 3.08567758149e18 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin UnitTemp_in_cgs: 1 # Kelvin
...@@ -11,17 +11,17 @@ TimeIntegration: ...@@ -11,17 +11,17 @@ TimeIntegration:
time_begin: 968 # The starting time of the simulation (in internal units). time_begin: 968 # The starting time of the simulation (in internal units).
time_end: 12000. # The end 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_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). dt_max: 10. # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics # Parameters governing the conserved quantities statistics
Statistics: Statistics:
delta_time: 1 # Time between statistics output delta_time: 24 # Time between statistics output
# Parameters governing the snapshots # Parameters governing the snapshots
Snapshots: Snapshots:
basename: Disc-Patch-dynamic # Common part of the name of output files basename: Disc-Patch-dynamic # Common part of the name of output files
time_first: 968. # Time of the first output (in internal units) time_first: 968. # Time of the first output (in internal units)
delta_time: 24. # Time difference between consecutive outputs (in internal units) delta_time: 96. # Time difference between outputs (in internal units)
# Parameters for the hydrodynamics scheme # Parameters for the hydrodynamics scheme
SPH: SPH:
...@@ -29,7 +29,7 @@ SPH: ...@@ -29,7 +29,7 @@ SPH:
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours. delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. 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_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). h_max: 60. # Maximal smoothing length allowed (in internal units).
# Parameters related to the initial conditions # Parameters related to the initial conditions
InitialConditions: InitialConditions:
...@@ -39,5 +39,7 @@ InitialConditions: ...@@ -39,5 +39,7 @@ InitialConditions:
DiscPatchPotential: DiscPatchPotential:
surface_density: 10. surface_density: 10.
scale_height: 100. scale_height: 100.
z_disc: 200. z_disc: 400.
z_trunc: 300.
z_max: 380.
timestep_mult: 0.03 timestep_mult: 0.03
examples/DiscPatch/HydroStatic/dynamic.pro deleted 100644 → 0
+ 0
− 139
View file @ 356a3348
;
; 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
examples/DiscPatch/HydroStatic/makeIC.py
+ 147
− 190
View file @ ad8aa1fb
...@@ -20,139 +20,147 @@ ...@@ -20,139 +20,147 @@
import h5py import h5py
import sys import sys
import numpy import numpy as np
import math import math
import random import random
import matplotlib.pyplot as plt
# Generates a disc-patch in hydrostatic equilibrium # Generates a disc-patch in hydrostatic equilibrium
# see Creasey, Theuns & Bower, 2013, for the equations: #
# disc parameters are: surface density sigma # See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
# scale height b #
# density: rho(z) = (sigma/2b) sech^2(z/b) #
# isothermal velocity dispersion = <v_z^2? = b pi G sigma # Disc parameters are: surface density -- sigma
# grad potential = 2 pi G sigma tanh(z/b) # scale height -- b
# potential = ln(cosh(z/b)) + const # gas adiabatic index -- gamma
# 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 # Problem parameters are: Ratio height/width of the box -- z_factor
# \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) # Size of the patch -- side_length
# usage: python makeIC.py 1000
# Parameters of the gas disc
surface_density = 10.
scale_height = 100.
gas_gamma = 5./3.
# Parameters of the problem
z_factor = 2
side_length = 400.
# File
fileName = "Disc-Patch.hdf5"
####################################################################
# physical constants in cgs # physical constants in cgs
NEWTON_GRAVITY_CGS = 6.672e-8 NEWTON_GRAVITY_CGS = 6.67408e-8
SOLAR_MASS_IN_CGS = 1.9885e33 SOLAR_MASS_IN_CGS = 1.9885e33
PARSEC_IN_CGS = 3.0856776e18 PARSEC_IN_CGS = 3.08567758149e18
PROTON_MASS_IN_CGS = 1.6726231e24 PROTON_MASS_IN_CGS = 1.672621898e-24
YEAR_IN_CGS = 3.154e+7 BOLTZMANN_IN_CGS = 1.38064852e-16
YEAR_IN_CGS = 3.15569252e7
# choice of units # choice of units
const_unit_length_in_cgs = (PARSEC_IN_CGS) unit_length_in_cgs = (PARSEC_IN_CGS)
const_unit_mass_in_cgs = (SOLAR_MASS_IN_CGS) unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
const_unit_velocity_in_cgs = (1e5) unit_velocity_in_cgs = (1e5)
unit_time_in_cgs = unit_length_in_cgs / unit_velocity_in_cgs
print "UnitMass_in_cgs: ", const_unit_mass_in_cgs
print "UnitLength_in_cgs: ", const_unit_length_in_cgs print "UnitMass_in_cgs: %.5e"%unit_mass_in_cgs
print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs print "UnitLength_in_cgs: %.5e"%unit_length_in_cgs
print "UnitVelocity_in_cgs: %.5e"%unit_velocity_in_cgs
print "UnitTime_in_cgs: %.5e"%unit_time_in_cgs
print ""
# Derived units
const_G = NEWTON_GRAVITY_CGS * unit_mass_in_cgs * unit_time_in_cgs**2 * unit_length_in_cgs**-3
const_mp = PROTON_MASS_IN_CGS * unit_mass_in_cgs**-1
const_kb = BOLTZMANN_IN_CGS * unit_mass_in_cgs**-1 * unit_length_in_cgs**-2 * unit_time_in_cgs**2
print "--- Some constants [internal units] ---"
print "G_Newton: %.5e"%const_G
print "m_proton: %.5e"%const_mp
print "k_boltzmann: %.5e"%const_kb
print ""
# derived quantities
temp = math.pi * const_G * surface_density * scale_height * const_mp / const_kb
u_therm = const_kb * temp / ((gas_gamma-1) * const_mp)
v_disp = math.sqrt(2 * u_therm)
soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma-1.)))
t_dyn = math.sqrt(scale_height / (const_G * surface_density))
t_cross = scale_height / soundspeed
print "--- Properties of the gas [internal units] ---"
print "Gas temperature: %.5e"%temp
print "Gas thermal_energy: %.5e"%u_therm
print "Dynamical time: %.5e"%t_dyn
print "Sound crossing time: %.5e"%t_cross
print "Gas sound speed: %.5e"%soundspeed
print "Gas 3D vel_disp: %.5e"%v_disp
print ""
# Problem properties
boxSize_x = side_length
boxSize_y = boxSize_x
boxSize_z = boxSize_x * z_factor
volume = boxSize_x * boxSize_y * boxSize_z
M_tot = boxSize_x * boxSize_y * surface_density * math.tanh(boxSize_z / (2. * scale_height))
density = M_tot / volume
entropy = (gas_gamma - 1.) * u_therm / density**(gas_gamma - 1.)
print "--- Problem properties [internal units] ---"
print "Box: [%.1f, %.1f, %.1f]"%(boxSize_x, boxSize_y, boxSize_z)
print "Volume: %.5e"%volume
print "Total mass: %.5e"%M_tot
print "Density: %.5e"%density
print "Entropy: %.5e"%entropy
print ""
####################################################################
# Read glass pre-ICs
infile = h5py.File('glassCube_32.hdf5', "r")
one_glass_pos = infile["/PartType0/Coordinates"][:,:]
one_glass_h = infile["/PartType0/SmoothingLength"][:]
# Rescale to the problem size
one_glass_pos *= boxSize_x
one_glass_h *= boxSize_x
#print min(one_glass_p[:,0]), max(one_glass_p[:,0])
#print min(one_glass_p[:,1]), max(one_glass_p[:,1])
#print min(one_glass_p[:,2]), max(one_glass_p[:,2])
# Now create enough copies to fill the volume in z
pos = np.copy(one_glass_pos)
h = np.copy(one_glass_h)
for i in range(1, z_factor):
one_glass_pos[:,2] += boxSize_x
pos = np.append(pos, one_glass_pos, axis=0)
h = np.append(h, one_glass_h, axis=0)
# parameters of potential #print min(pos[:,0]), max(pos[:,0])
surface_density = 100. # surface density of all mass, which generates the gravitational potential #print min(pos[:,1]), max(pos[:,1])
scale_height = 100. #print min(pos[:,2]), max(pos[:,2])
gamma = 5./3.
fgas = 0.1 # gas fraction
# 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
# Compute further properties of ICs
numPart = np.size(h)
mass = M_tot / numPart
# Parameters print "--- Particle properties [internal units] ---"
periodic= 1 # 1 For periodic box print "Number part.: ", numPart
boxSize = 400. # [kpc] print "Part. mass: %.5e"%mass
Radius = 100. # maximum radius of particles [kpc] print ""
G = const_G
# File # Create additional arrays
fileName = "Disc-Patch.hdf5" u = np.ones(numPart) * u_therm
mass = np.ones(numPart) * mass
vel = np.zeros((numPart, 3))
ids = 1 + np.linspace(0, numPart, numPart, endpoint=False)
#---------------------------------------------------
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 = fgas * 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 disc to middle of box
pos[:,:] += boxSize/2
# create numPart dm particles
numPart = 0
####################################################################
# Create and write output file # Create and write output file
#File #File
...@@ -160,97 +168,46 @@ file = h5py.File(fileName, 'w') ...@@ -160,97 +168,46 @@ file = h5py.File(fileName, 'w')
#Units #Units
grp = file.create_group("/Units") grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs grp.attrs["Unit length in cgs (U_L)"] = unit_length_in_cgs
grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs grp.attrs["Unit mass in cgs (U_M)"] = 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 time in cgs (U_t)"] = unit_time_in_cgs
grp.attrs["Unit current in cgs (U_I)"] = 1. grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1. grp.attrs["Unit temperature in cgs (U_T)"] = 1.
# Header # Header
grp = file.create_group("/Header") grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize grp.attrs["BoxSize"] = [boxSize_x, boxSize_y, boxSize_z]
grp.attrs["NumPart_Total"] = [numGas, numPart, 0, 0, 0, 0] grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 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["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0 grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1 grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = [entropy_flag] grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
grp.attrs["Dimension"] = 3 grp.attrs["Dimension"] = 3
#Runtime parameters #Runtime parameters
grp = file.create_group("/RuntimePars") grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic grp.attrs["PeriodicBoundariesOn"] = 1
# write gas particles # write gas particles
grp0 = file.create_group("/PartType0") grp0 = file.create_group("/PartType0")
ds = grp0.create_dataset('Coordinates', (numGas, 3), 'f') ds = grp0.create_dataset('Coordinates', (numPart, 3), 'f', data=pos)
ds[()] = pos ds = grp0.create_dataset('Velocities', (numPart, 3), 'f')
ds = grp0.create_dataset('Masses', (numPart,), 'f', data=mass)
ds = grp0.create_dataset('SmoothingLength', (numPart,), 'f', data=h)
ds = grp0.create_dataset('InternalEnergy', (numPart,), 'f', data=u)
ds = grp0.create_dataset('ParticleIDs', (numPart, ), 'L', data=ids)
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)
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() print "--- Runtime parameters (YAML file): ---"
print "DiscPatchPotential:surface_density: ", surface_density
print "DiscPatchPotential:scale_height: ", scale_height
print "DiscPatchPotential:z_disc: ", boxSize_z / 2.
print ""
sys.exit() print "--- Constant parameters: ---"
print "const_isothermal_internal_energy: %ef"%u_therm
examples/DiscPatch/HydroStatic/plot.py → examples/DiscPatch/HydroStatic/plotSolution.py
+ 23
− 6
View file @ ad8aa1fb
...@@ -20,7 +20,7 @@ ...@@ -20,7 +20,7 @@
## ##
# This script plots the Disc-Patch_*.hdf5 snapshots. # This script plots the Disc-Patch_*.hdf5 snapshots.
# It takes two (optional) parameters: the counter value of the first and last # 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 import numpy as np
...@@ -34,12 +34,14 @@ import sys ...@@ -34,12 +34,14 @@ import sys
# Parameters # Parameters
surface_density = 10. surface_density = 10.
scale_height = 100. scale_height = 100.
z_disc = 200. z_disc = 400.
utherm = 20.2615290634 z_trunc = 300.
z_max = 350.
utherm = 20.2678457288
gamma = 5. / 3. gamma = 5. / 3.
start = 0 start = 0
stop = 81 stop = 21
if len(sys.argv) > 1: if len(sys.argv) > 1:
start = int(sys.argv[1]) start = int(sys.argv[1])
if len(sys.argv) > 2: if len(sys.argv) > 2:
...@@ -78,22 +80,37 @@ for f in sorted(glob.glob("Disc-Patch_*.hdf5")): ...@@ -78,22 +80,37 @@ for f in sorted(glob.glob("Disc-Patch_*.hdf5")):
print "processing", f, "..." print "processing", f, "..."
zrange = np.linspace(0., 400., 1000) zrange = np.linspace(0., 2. * z_disc, 1000)
time, z, rho, P, v = get_data(f) time, z, rho, P, v = get_data(f)
fig, ax = pl.subplots(3, 1, sharex = True) fig, ax = pl.subplots(3, 1, sharex = True)
ax[0].plot(z, rho, "r.") ax[0].plot(z, rho, "r.")
ax[0].plot(zrange, get_analytic_density(zrange), "k-") ax[0].plot(zrange, get_analytic_density(zrange), "k-")
ax[0].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
ax[0].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
ax[0].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
ax[0].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
ax[0].set_ylim(0., 1.2 * get_analytic_density(z_disc))
ax[0].set_ylabel("density") ax[0].set_ylabel("density")
ax[1].plot(z, v, "r.") ax[1].plot(z, v, "r.")
ax[1].plot(zrange, np.zeros(len(zrange)), "k-") ax[1].plot(zrange, np.zeros(len(zrange)), "k-")
ax[1].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
ax[1].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
ax[1].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
ax[1].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
ax[1].set_ylim(-0.5, 10.)
ax[1].set_ylabel("velocity norm") ax[1].set_ylabel("velocity norm")
ax[2].plot(z, P, "r.") ax[2].plot(z, P, "r.")
ax[2].plot(zrange, get_analytic_pressure(zrange), "k-") ax[2].plot(zrange, get_analytic_pressure(zrange), "k-")
ax[2].set_xlim(0., 400.) ax[2].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
ax[2].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
ax[2].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
ax[2].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
ax[2].set_xlim(0., 2. * z_disc)
ax[2].set_ylim(0., 1.2 * get_analytic_pressure(z_disc))
ax[2].set_xlabel("z") ax[2].set_xlabel("z")
ax[2].set_ylabel("pressure") ax[2].set_ylabel("pressure")
......
examples/DiscPatch/HydroStatic/test.pro deleted 100644 → 0
+ 0
− 142
View file @ 356a3348
;
; 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
examples/parameter_example.yml
+ 3
− 1
View file @ ad8aa1fb
...@@ -103,7 +103,9 @@ IsothermalPotential: ...@@ -103,7 +103,9 @@ IsothermalPotential:
DiscPatchPotential: DiscPatchPotential:
surface_density: 10. # Surface density of the disc (internal units) surface_density: 10. # Surface density of the disc (internal units)
scale_height: 100. # Scale height of the disc (internal units) scale_height: 100. # Scale height of the disc (internal units)
z_disc: 200. # Position of the disc along the z-axis (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.
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition 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) growth_time: 5. # (Optional) Time for the disc to grow to its final size (multiple of the dynamical time)
......
src/const.h
+ 1
− 1
View file @ ad8aa1fb
...@@ -37,7 +37,7 @@ ...@@ -37,7 +37,7 @@
#define const_max_u_change 0.1f #define const_max_u_change 0.1f
/* Thermal energy per unit mass used as a constant for the isothermal EoS */ /* Thermal energy per unit mass used as a constant for the isothermal EoS */
#define const_isothermal_internal_energy 20.2615290634f #define const_isothermal_internal_energy 20.2678457288f
/* Type of gradients to use (GIZMO_SPH only) */ /* Type of gradients to use (GIZMO_SPH only) */
/* If no option is chosen, no gradients are used (first order scheme) */ /* If no option is chosen, no gradients are used (first order scheme) */
......
src/potential/disc_patch/potential.h
+ 162
− 49
View file @ ad8aa1fb
...@@ -39,34 +39,63 @@ ...@@ -39,34 +39,63 @@
/** /**
* @brief External Potential Properties - Disc patch case * @brief External Potential Properties - Disc patch case
* *
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948 * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948.
*
* We truncate the accelerations beyond z_trunc using a 1-cos(z) function
* that smoothly brings the accelerations to 0 at z_max.
*/ */
struct external_potential { struct external_potential {
/*! Surface density of the disc */ /*! Surface density of the disc (sigma) */
double surface_density; float surface_density;
/*! Disc scale-height (b) */
float scale_height;
/*! Disc scale-height */ /*! Inverse of disc scale-height (1/b) */
double scale_height; float scale_height_inv;
/*! Position of the disc along the z-axis */ /*! Position of the disc along the z-axis */
double z_disc; float z_disc;
/*! Position above which the accelerations get truncated */
float z_trunc;
/*! Position above which the accelerations are zero */
float z_max;
/*! The truncated transition regime */
float z_trans;
/*! Inverse of the truncated transition regime */
float z_trans_inv;
/*! Dynamical time of the system */ /*! Dynamical time of the system */
double dynamical_time; float dynamical_time;
/*! Time over which to grow the disk in units of the dynamical time */ /*! Time over which to grow the disk */
double growth_time; float growth_time;
/*! Inverse of the growth time */
float growth_time_inv;
/*! Time-step condition pre-factor */ /*! Time-step condition pre-factor */
double timestep_mult; float timestep_mult;
/*! Constant pre-factor (2 pi G sigma) */
float norm;
/*! Constant pre-factor (2 pi sigma)*/
float norm_over_G;
}; };
/** /**
* @brief Computes the time-step from the acceleration due to a hydrostatic * @brief Computes the time-step from the acceleration due to a hydrostatic
* disc. * disc.
* *
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948 * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948,
* equations 17 and 20.
* We do not use the truncated potential here.
* *
* @param time The current time. * @param time The current time.
* @param potential The properties of the potential. * @param potential The properties of the potential.
...@@ -80,39 +109,41 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( ...@@ -80,39 +109,41 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
/* initilize time step to disc dynamical time */ /* initilize time step to disc dynamical time */
const float dt_dyn = potential->dynamical_time; const float dt_dyn = potential->dynamical_time;
float dt = dt_dyn; const float b = potential->scale_height;
const float b_inv = potential->scale_height_inv;
const float norm = potential->norm;
/* absolute value of height above disc */ /* absolute value of height above disc */
const float dz = fabsf(g->x[2] - potential->z_disc); const float dz = fabsf(g->x[2] - potential->z_disc);
/* vertical acceleration */ /* vertical acceleration */
const float z_accel = 2.f * M_PI * phys_const->const_newton_G * const float z_accel = norm * tanhf(dz * b_inv);
potential->surface_density *
tanhf(dz / potential->scale_height); float dt = dt_dyn;
/* demand that dt * velocity < fraction of scale height of disc */ /* demand that dt * velocity < fraction of scale height of disc */
float dt1 = FLT_MAX;
if (g->v_full[2] != 0.f) { if (g->v_full[2] != 0.f) {
dt1 = potential->scale_height / fabsf(g->v_full[2]);
if (dt1 < dt) dt = dt1; const float dt1 = b / fabsf(g->v_full[2]);
dt = min(dt1, dt);
} }
/* demand that dt^2 * acceleration < fraction of scale height of disc */ /* demand that dt^2 * acceleration < fraction of scale height of disc */
float dt2 = FLT_MAX;
if (z_accel != 0.f) { if (z_accel != 0.f) {
dt2 = potential->scale_height / fabsf(z_accel);
const float dt2 = b / fabsf(z_accel);
if (dt2 < dt * dt) dt = sqrtf(dt2); if (dt2 < dt * dt) dt = sqrtf(dt2);
} }
/* demand that dt^3 * jerk < fraction of scale height of disc */ /* demand that dt^3 * jerk < fraction of scale height of disc */
float dt3 = FLT_MAX;
if (g->v_full[2] != 0.f) { if (g->v_full[2] != 0.f) {
const float cosh_dz_inv = 1.f / coshf(dz * b_inv);
const float cosh_dz_inv2 = cosh_dz_inv * cosh_dz_inv;
const float dz_accel_over_dt = const float dz_accel_over_dt =
2.f * M_PI * phys_const->const_newton_G * potential->surface_density / norm * cosh_dz_inv2 * b_inv * fabsf(g->v_full[2]);
potential->scale_height / coshf(dz / potential->scale_height) /
coshf(dz / potential->scale_height) * fabsf(g->v_full[2]);
dt3 = potential->scale_height / fabsf(dz_accel_over_dt); const float dt3 = b / fabsf(dz_accel_over_dt);
if (dt3 < dt * dt * dt) dt = cbrtf(dt3); if (dt3 < dt * dt * dt) dt = cbrtf(dt3);
} }
...@@ -125,6 +156,8 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( ...@@ -125,6 +156,8 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
* *
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948, * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948,
* equation 17. * equation 17.
* We truncate the accelerations beyond z_trunc using a 1-cos(z) function
* that smoothly brings the accelerations to 0 at z_max.
* *
* @param time The current time in internal units. * @param time The current time in internal units.
* @param potential The properties of the potential. * @param potential The properties of the potential.
...@@ -136,19 +169,38 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( ...@@ -136,19 +169,38 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
const struct phys_const* restrict phys_const, struct gpart* restrict g) { const struct phys_const* restrict phys_const, struct gpart* restrict g) {
const float dz = g->x[2] - potential->z_disc; const float dz = g->x[2] - potential->z_disc;
const float t_dyn = potential->dynamical_time; const float abs_dz = fabsf(dz);
const float t_growth = potential->growth_time;
float reduction_factor = 1.f; const float t_growth_inv = potential->growth_time_inv;
if (time < potential->growth_time * t_dyn) const float b_inv = potential->scale_height_inv;
reduction_factor = time / (potential->growth_time * t_dyn); const float z_trunc = potential->z_trunc;
const float z_max = potential->z_max;
/* Accelerations. Note that they are multiplied by G later on */ const float z_trans_inv = potential->z_trans_inv;
const float z_accel = reduction_factor * 2.f * M_PI * const float norm_over_G = potential->norm_over_G;
potential->surface_density *
tanhf(fabsf(dz) / potential->scale_height); /* Are we still growing the disc ? */
const float reduction_factor = time < t_growth ? time * t_growth_inv : 1.f;
/* Truncated or not ? */
float a_z;
if (abs_dz < z_trunc) {
/* Acc. 2 pi sigma tanh(z/b) */
a_z = reduction_factor * norm_over_G * tanhf(abs_dz * b_inv);
} else if (abs_dz < z_max) {
/* Acc. 2 pi sigma tanh(z/b) [1/2 + 1/2cos((z-zmax)/(pi z_trans))] */
a_z = reduction_factor * norm_over_G * tanhf(abs_dz * b_inv) *
(0.5f + 0.5f * cosf((float)(M_PI) * (abs_dz - z_trunc) * z_trans_inv));
} else {
/* Acc. 0 */
a_z = 0.f;
}
if (dz > 0) g->a_grav[2] -= z_accel; /* Get the correct sign. Recall G is multipiled in later on */
if (dz < 0) g->a_grav[2] += z_accel; if (dz > 0) g->a_grav[2] -= a_z;
if (dz < 0) g->a_grav[2] += a_z;
} }
/** /**
...@@ -156,7 +208,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( ...@@ -156,7 +208,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
* disc patch potential. * disc patch potential.
* *
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948, * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948,
* equation 24. * equation 22.
* We truncate the accelerations beyond z_trunc using a 1-cos(z) function
* that smoothly brings the accelerations to 0 at z_max.
* *
* @param time The current time. * @param time The current time.
* @param potential The #external_potential used in the run. * @param potential The #external_potential used in the run.
...@@ -169,16 +223,35 @@ external_gravity_get_potential_energy( ...@@ -169,16 +223,35 @@ external_gravity_get_potential_energy(
const struct phys_const* const phys_const, const struct gpart* gp) { const struct phys_const* const phys_const, const struct gpart* gp) {
const float dz = gp->x[2] - potential->z_disc; const float dz = gp->x[2] - potential->z_disc;
const float t_dyn = potential->dynamical_time; const float abs_dz = fabsf(dz);
const float t_growth = potential->growth_time;
const float t_growth_inv = potential->growth_time_inv;
const float b = potential->scale_height;
const float b_inv = potential->scale_height_inv;
const float norm = potential->norm;
const float z_trunc = potential->z_trunc;
const float z_max = potential->z_max;
/* Are we still growing the disc ? */
const float reduction_factor = time < t_growth ? time * t_growth_inv : 1.f;
/* Truncated or not ? */
float pot;
if (abs_dz < z_trunc) {
/* Potential (2 pi G sigma b ln(cosh(z/b)) */
pot = b * logf(coshf(dz * b_inv));
} else if (abs_dz < z_max) {
/* Potential. At z>>b, phi(z) = norm * z / b */
pot = 0.f;
float reduction_factor = 1.f; } else {
if (time < potential->growth_time * t_dyn)
reduction_factor = time / (potential->growth_time * t_dyn);
/* Accelerations. Note that they are multiplied by G later on */ pot = 0.f;
return reduction_factor * 2.f * M_PI * phys_const->const_newton_G * }
potential->surface_density * potential->scale_height *
logf(coshf(dz / potential->scale_height)); return pot * reduction_factor * norm;
} }
/** /**
...@@ -204,13 +277,47 @@ static INLINE void potential_init_backend( ...@@ -204,13 +277,47 @@ static INLINE void potential_init_backend(
parameter_file, "DiscPatchPotential:scale_height"); parameter_file, "DiscPatchPotential:scale_height");
potential->z_disc = potential->z_disc =
parser_get_param_double(parameter_file, "DiscPatchPotential:z_disc"); parser_get_param_double(parameter_file, "DiscPatchPotential:z_disc");
potential->z_trunc = parser_get_opt_param_double(
parameter_file, "DiscPatchPotential:z_trunc", FLT_MAX);
potential->z_max = parser_get_opt_param_double(
parameter_file, "DiscPatchPotential:z_max", FLT_MAX);
potential->z_disc =
parser_get_param_double(parameter_file, "DiscPatchPotential:z_disc");
potential->timestep_mult = parser_get_param_double( potential->timestep_mult = parser_get_param_double(
parameter_file, "DiscPatchPotential:timestep_mult"); parameter_file, "DiscPatchPotential:timestep_mult");
potential->growth_time = parser_get_opt_param_double( potential->growth_time = parser_get_opt_param_double(
parameter_file, "DiscPatchPotential:growth_time", 0.); parameter_file, "DiscPatchPotential:growth_time", 0.);
/* Compute the dynamical time */
potential->dynamical_time = potential->dynamical_time =
sqrt(potential->scale_height / sqrt(potential->scale_height /
(phys_const->const_newton_G * potential->surface_density)); (phys_const->const_newton_G * potential->surface_density));
/* Convert the growth time multiplier to physical time */
potential->growth_time *= potential->dynamical_time;
/* Some cross-checks */
if (potential->z_trunc > potential->z_max)
error("Potential truncation z larger than maximal z");
if (potential->z_trunc < potential->scale_height)
error("Potential truncation z smaller than scale height");
/* Compute derived quantities */
potential->scale_height_inv = 1. / potential->scale_height;
potential->norm =
2. * M_PI * phys_const->const_newton_G * potential->surface_density;
potential->norm_over_G = 2 * M_PI * potential->surface_density;
potential->z_trans = potential->z_max - potential->z_trunc;
if (potential->z_trans != 0.f)
potential->z_trans_inv = 1. / potential->z_trans;
else
potential->z_trans_inv = FLT_MAX;
if (potential->growth_time != 0.)
potential->growth_time_inv = 1. / potential->growth_time;
else
potential->growth_time_inv = FLT_MAX;
} }
/** /**
...@@ -222,13 +329,19 @@ static INLINE void potential_print_backend( ...@@ -222,13 +329,19 @@ static INLINE void potential_print_backend(
const struct external_potential* potential) { const struct external_potential* potential) {
message( message(
"External potential is 'Disk-patch' with properties surface_density = %e " "External potential is 'Disk-patch' with Sigma=%f, z_disc=%f, b=%f and "
"disc height= %e scale height = %e timestep multiplier = %e.", "dt_mult=%f.",
potential->surface_density, potential->z_disc, potential->scale_height, potential->surface_density, potential->z_disc, potential->scale_height,
potential->timestep_mult); potential->timestep_mult);
if (potential->z_max < FLT_MAX)
message("Potential will be truncated at z_trunc=%f and zeroed at z_max=%f",
potential->z_trunc, potential->z_max);
if (potential->growth_time > 0.) if (potential->growth_time > 0.)
message("Disc will grow for %f dynamical times.", potential->growth_time); message("Disc will grow for %f [time_units]. (%f dynamical time)",
potential->growth_time,
potential->growth_time / potential->dynamical_time);
} }
#endif /* SWIFT_DISC_PATCH_H */ #endif /* SWIFT_DISC_PATCH_H */
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment