diff --git a/examples/DiscPatch/HydroStatic/disc-patch-icc.yml b/examples/DiscPatch/HydroStatic/disc-patch-icc.yml index 6a27016b8a3f484b7c1c9b74594073d5f28efe90..6a5735f06ed4c0ea9884fabb793c9d2d84bd8105 100644 --- a/examples/DiscPatch/HydroStatic/disc-patch-icc.yml +++ b/examples/DiscPatch/HydroStatic/disc-patch-icc.yml @@ -1,8 +1,8 @@ # 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 + UnitMass_in_cgs: 1.9885e33 # Grams + UnitLength_in_cgs: 3.08567758149e18 # Centimeters + UnitVelocity_in_cgs: 1e5 # Centimeters per second UnitCurrent_in_cgs: 1 # Amperes UnitTemp_in_cgs: 1 # Kelvin @@ -11,17 +11,17 @@ 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). + dt_max: 10. # 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 + delta_time: 12. # Time between statistics output # Parameters governing the snapshots Snapshots: - basename: Disc-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) + basename: Disc-Patch # Common part of the name of output files + time_first: 0. # Time of the first output (in internal units) + delta_time: 48. # Time difference between outputs (in internal units) # Parameters for the hydrodynamics scheme SPH: @@ -29,7 +29,7 @@ SPH: 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). + h_max: 60. # Maximal smoothing length allowed (in internal units). # Parameters related to the initial conditions InitialConditions: @@ -39,6 +39,6 @@ InitialConditions: DiscPatchPotential: surface_density: 10. scale_height: 100. - z_disc: 200. + z_disc: 400. timestep_mult: 0.03 growth_time: 5. diff --git a/examples/DiscPatch/HydroStatic/disc-patch.yml b/examples/DiscPatch/HydroStatic/disc-patch.yml index 8bd67c5b08de82bb6a3d47ccf3419f85e3e5c6b1..9454cadf38afb3ea6bfa6c978b6894261cfdf021 100644 --- a/examples/DiscPatch/HydroStatic/disc-patch.yml +++ b/examples/DiscPatch/HydroStatic/disc-patch.yml @@ -1,8 +1,8 @@ # 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 + UnitMass_in_cgs: 1.9885e33 # Grams + UnitLength_in_cgs: 3.08567758149e18 # Centimeters + UnitVelocity_in_cgs: 1e5 # Centimeters per second UnitCurrent_in_cgs: 1 # Amperes UnitTemp_in_cgs: 1 # Kelvin @@ -11,17 +11,17 @@ TimeIntegration: time_begin: 968 # The starting 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_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 Statistics: - delta_time: 1 # Time between statistics output + delta_time: 24 # Time between statistics output # Parameters governing the snapshots Snapshots: - basename: Disc-Patch-dynamic # Common part of the name of output files - time_first: 968. # Time of the first output (in internal units) - delta_time: 24. # Time difference between consecutive outputs (in internal units) + basename: Disc-Patch-dynamic # Common part of the name of output files + time_first: 968. # Time of the first output (in internal units) + delta_time: 96. # Time difference between outputs (in internal units) # Parameters for the hydrodynamics scheme SPH: @@ -29,7 +29,7 @@ SPH: 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). + h_max: 60. # Maximal smoothing length allowed (in internal units). # Parameters related to the initial conditions InitialConditions: @@ -39,5 +39,5 @@ InitialConditions: DiscPatchPotential: surface_density: 10. scale_height: 100. - z_disc: 200. + z_disc: 400. timestep_mult: 0.03 diff --git a/examples/DiscPatch/HydroStatic/makeIC.py b/examples/DiscPatch/HydroStatic/makeIC.py index 6ba1ccd06fed84ca728aadaa5922dbba536b6881..423dc8a44714a6e710b0c037f8a33d2120fbdcc8 100644 --- a/examples/DiscPatch/HydroStatic/makeIC.py +++ b/examples/DiscPatch/HydroStatic/makeIC.py @@ -20,139 +20,147 @@ import h5py import sys -import numpy +import numpy as np import math import random -import matplotlib.pyplot as plt # Generates a disc-patch in hydrostatic equilibrium -# 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 -# grad potential = 2 pi G sigma tanh(z/b) -# potential = ln(cosh(z/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) -# usage: python makeIC.py 1000 +# +# See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948 +# +# +# Disc parameters are: surface density -- sigma +# scale height -- b +# gas adiabatic index -- gamma +# +# Problem parameters are: Ratio height/width of the box -- z_factor +# Size of the patch -- side_length + +# 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 -NEWTON_GRAVITY_CGS = 6.672e-8 +NEWTON_GRAVITY_CGS = 6.67408e-8 SOLAR_MASS_IN_CGS = 1.9885e33 -PARSEC_IN_CGS = 3.0856776e18 -PROTON_MASS_IN_CGS = 1.6726231e24 -YEAR_IN_CGS = 3.154e+7 +PARSEC_IN_CGS = 3.08567758149e18 +PROTON_MASS_IN_CGS = 1.672621898e-24 +BOLTZMANN_IN_CGS = 1.38064852e-16 +YEAR_IN_CGS = 3.15569252e7 # 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) - -print "UnitMass_in_cgs: ", const_unit_mass_in_cgs -print "UnitLength_in_cgs: ", const_unit_length_in_cgs -print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs +unit_length_in_cgs = (PARSEC_IN_CGS) +unit_mass_in_cgs = (SOLAR_MASS_IN_CGS) +unit_velocity_in_cgs = (1e5) +unit_time_in_cgs = unit_length_in_cgs / unit_velocity_in_cgs + +print "UnitMass_in_cgs: %.5e"%unit_mass_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 -surface_density = 100. # surface density of all mass, which generates the gravitational potential -scale_height = 100. -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 +#print min(pos[:,0]), max(pos[:,0]) +#print min(pos[:,1]), max(pos[:,1]) +#print min(pos[:,2]), max(pos[:,2]) +# Compute further properties of ICs +numPart = np.size(h) +mass = M_tot / numPart -# Parameters -periodic= 1 # 1 For periodic box -boxSize = 400. # [kpc] -Radius = 100. # maximum radius of particles [kpc] -G = const_G +print "--- Particle properties [internal units] ---" +print "Number part.: ", numPart +print "Part. mass: %.5e"%mass +print "" -# File -fileName = "Disc-Patch.hdf5" +# Create additional arrays +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 #File @@ -160,97 +168,46 @@ file = h5py.File(fileName, 'w') #Units grp = file.create_group("/Units") -grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs -grp.attrs["Unit mass in cgs (U_M)"] = const_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 length in cgs (U_L)"] = unit_length_in_cgs +grp.attrs["Unit mass in cgs (U_M)"] = unit_mass_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 temperature in cgs (U_T)"] = 1. # Header grp = file.create_group("/Header") -grp.attrs["BoxSize"] = boxSize -grp.attrs["NumPart_Total"] = [numGas, numPart, 0, 0, 0, 0] +grp.attrs["BoxSize"] = [boxSize_x, boxSize_y, boxSize_z] +grp.attrs["NumPart_Total"] = [numPart, 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["NumFilesPerSnapshot"] = 1 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 #Runtime parameters grp = file.create_group("/RuntimePars") -grp.attrs["PeriodicBoundariesOn"] = periodic - +grp.attrs["PeriodicBoundariesOn"] = 1 # write gas particles grp0 = file.create_group("/PartType0") -ds = grp0.create_dataset('Coordinates', (numGas, 3), 'f') -ds[()] = pos +ds = grp0.create_dataset('Coordinates', (numPart, 3), 'f', data=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 diff --git a/src/const.h b/src/const.h index 141eb48acc633542aa98655caa8debdd2dbce530..f5232de0b2b4aba6701c07f1e5f91d5be3d5b6cc 100644 --- a/src/const.h +++ b/src/const.h @@ -37,7 +37,7 @@ #define const_max_u_change 0.1f /* 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) */ /* If no option is chosen, no gradients are used (first order scheme) */ diff --git a/src/potential/disc_patch/potential.h b/src/potential/disc_patch/potential.h index 8fa40ecd4e6503cde8be00db8c6fb8a70c84ebdf..5d729825b647962a8cb9c94ab1d7bafc1a9190d7 100644 --- a/src/potential/disc_patch/potential.h +++ b/src/potential/disc_patch/potential.h @@ -66,7 +66,8 @@ struct external_potential { * @brief Computes the time-step from the acceleration due to a hydrostatic * 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. * * @param time The current time. * @param potential The properties of the potential. @@ -156,7 +157,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( * disc patch potential. * * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948, - * equation 24. + * equation 22. * * @param time The current time. * @param potential The #external_potential used in the run.