Updated the hydrodstatic disc patch example to use a longer slab in z and the...

Updated the hydrodstatic disc patch example to use a longer slab in z and the correct internal energy.

examples/DiscPatch/HydroStatic/disc-patch-icc.yml

 # 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.

examples/DiscPatch/HydroStatic/disc-patch.yml

 # 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

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

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) */

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.