diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml index f77b4632c67115559ab2d45083eb7778d21903d8..20b5bb3aa613d553d8c401e968d8ebfc0572e610 100644 --- a/examples/ExternalPointMass/externalPointMass.yml +++ b/examples/ExternalPointMass/externalPointMass.yml @@ -31,7 +31,7 @@ SPH: # Parameters related to the initial conditions InitialConditions: - file_name: Sphere.hdf5 # The file to read + file_name: PointMass.hdf5 # The file to read shift_x: 50. # A shift to apply to all particles read from the ICs (in internal units). shift_y: 50. shift_z: 50. diff --git a/examples/ExternalPointMass/makeIC.py b/examples/ExternalPointMass/makeIC.py index a324220e183baed4ac5aaaf2d71cf5b6ea15c1ce..8bda95681a6e9b71d54097f5cba8efd6ad2bd537 100644 --- a/examples/ExternalPointMass/makeIC.py +++ b/examples/ExternalPointMass/makeIC.py @@ -24,7 +24,7 @@ import numpy import math import random -# Generates a random distriution of particles, for motion in an external potnetial centred at (0,0,0) +# Generates a random distriution of particles, for motion in an external potential centred at (0,0,0) # physical constants in cgs NEWTON_GRAVITY_CGS = 6.67408e-8 @@ -39,34 +39,28 @@ 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 +print "UnitTime_in_cgs: ", const_unit_length_in_cgs / const_unit_velocity_in_cgs # 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 +print '---------------------' +print 'G in internal units: ', const_G # Parameters periodic= 1 # 1 For periodic box boxSize = 100. # Radius = boxSize / 4. # maximum radius of particles -G = const_G Mass = 1e10 +print "Mass at the centre: ", Mass -N = int(sys.argv[1]) # Number of particles -L = N**(1./3.) +numPart = int(sys.argv[1]) # Number of particles +mass = 1. -# these are not used but necessary for I/O -rho = 2. # Density -P = 1. # Pressure -gamma = 5./3. # Gas adiabatic index -fileName = "Sphere.hdf5" +fileName = "PointMass.hdf5" -#--------------------------------------------------- -numPart = N -mass = 1 -internalEnergy = P / ((gamma - 1.)*rho) #-------------------------------------------------- @@ -98,24 +92,24 @@ grp.attrs["Unit current in cgs (U_I)"] = 1. grp.attrs["Unit temperature in cgs (U_T)"] = 1. #Particle group -#grp0 = file.create_group("/PartType0") grp1 = file.create_group("/PartType1") + #generate particle positions -radius = Radius * (numpy.random.rand(N))**(1./3.) -ctheta = -1. + 2 * numpy.random.rand(N) +radius = Radius * (numpy.random.rand(numPart))**(1./3.) +ctheta = -1. + 2 * numpy.random.rand(numPart) stheta = numpy.sqrt(1.-ctheta**2) -phi = 2 * math.pi * numpy.random.rand(N) +phi = 2 * math.pi * numpy.random.rand(numPart) r = numpy.zeros((numPart, 3)) -# r[:,0] = radius * stheta * numpy.cos(phi) -# r[:,1] = radius * stheta * numpy.sin(phi) -# r[:,2] = radius * ctheta r[:,0] = radius -# -speed = numpy.sqrt(G * Mass / radius) + +#generate particle velocities +speed = numpy.sqrt(const_G * Mass / radius) v = numpy.zeros((numPart, 3)) omega = speed / radius period = 2.*math.pi/omega -print 'period = minimum = ',min(period), ' maximum = ',max(period) +print '---------------------' +print 'Period: minimum = ',min(period) +print 'Period: maximum = ',max(period) v[:,0] = -omega * r[:,1] v[:,1] = omega * r[:,0] @@ -129,17 +123,6 @@ ds = grp1.create_dataset('Masses', (numPart,), 'f') ds[()] = m m = numpy.zeros(1) -h = numpy.full((numPart, ), 1.1255 * boxSize / L) -ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f') -ds[()] = h -h = numpy.zeros(1) - -u = numpy.full((numPart, ), internalEnergy) -ds = grp1.create_dataset('InternalEnergy', (numPart,), 'f') -ds[()] = u -u = numpy.zeros(1) - - ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False) ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L') ds[()] = ids diff --git a/examples/ExternalPointMass/run.sh b/examples/ExternalPointMass/run.sh index 44daaaaf5b92c2923527daf256a61fbb4403b73f..e074c384c4e002a161c7d8258e9068663204099f 100755 --- a/examples/ExternalPointMass/run.sh +++ b/examples/ExternalPointMass/run.sh @@ -1,7 +1,7 @@ #!/bin/bash # Generate the initial conditions if they are not present. -if [ ! -e Sphere.hdf5 ] +if [ ! -e PointMass.hdf5 ] then echo "Generating initial conditions for the point mass potential box example..." python makeIC.py 10000