From cc2dfa53ae9a3ce4fd7e20f0ab54ee484a6e0f31 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Tue, 7 Feb 2017 12:30:28 +0000
Subject: [PATCH] Removed unnecessary lines in script to generate ICs for
external point mass example.
---
.../ExternalPointMass/externalPointMass.yml | 2 +-
examples/ExternalPointMass/makeIC.py | 53 +++++++------------
examples/ExternalPointMass/run.sh | 2 +-
3 files changed, 20 insertions(+), 37 deletions(-)
diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml
index f77b4632c6..20b5bb3aa6 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 a324220e18..8bda95681a 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 44daaaaf5b..e074c384c4 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
--
GitLab