Removed unwanted changes to the Isothermal potential test-case

examples/IsothermalPotential/GravityOnly/isothermal.yml

@@ -9,7 +9,7 @@ InternalUnitSystem:
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   10.    # The end time of the simulation (in internal units).
+  time_end:   8.    # The end time of the simulation (in internal units).
   dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-1  # The maximal time-step size of the simulation (in internal units).
 
@@ -21,7 +21,7 @@ Statistics:
 Snapshots:
   basename:            Isothermal # Common part of the name of output files
   time_first:          0.         # Time of the first output (in internal units)
-  delta_time:          0.1       # Time difference between consecutive outputs (in internal units)
+  delta_time:          0.02       # Time difference between consecutive outputs (in internal units)
 
 # Parameters for the hydrodynamics scheme
 SPH:
@@ -33,15 +33,15 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  Isothermal.hdf5       # The file to read
-  shift_x:    200.                  # A shift to apply to all particles read from the ICs (in internal units).
-  shift_y:    200.
-  shift_z:    200.
+  shift_x:    100.                  # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    100.
+  shift_z:    100.
  
 # External potential parameters
 IsothermalPotential:
-  position_x:      200.     # location of centre of isothermal potential in internal units
-  position_y:      200.
-  position_z:      200.	
+  position_x:      100.     # location of centre of isothermal potential in internal units
+  position_y:      100.
+  position_z:      100.	
   vrot:            200.     # rotation speed of isothermal potential in internal units
   timestep_mult:   0.03     # controls time step
 

examples/IsothermalPotential/GravityOnly/makeIC.py

@@ -25,8 +25,8 @@ import math
 import random
 
 # Generates N particles in a spherical distribution centred on [0,0,0], to be moved in an isothermal potential
-# usage: python makeIC.py 1000 0 1: generate 1000 particles on circular orbits
-#        python makeIC.py 1000 1 1: generate 1000 particles with Lz/L uniform in [0,1]
+# usage: python makeIC.py 1000 0 : generate 1000 particles on circular orbits
+#        python makeIC.py 1000 1 : generate 1000 particles with Lz/L uniform in [0,1]
 # all particles move in the xy plane, and start at y=0
 
 # physical constants in cgs
@@ -64,7 +64,6 @@ G       = const_G
 
 N       = int(sys.argv[1])  # Number of particles
 icirc   = int(sys.argv[2])  # if = 0, all particles are on circular orbits, if = 1, Lz/Lcirc uniform in ]0,1[
-three_d = int(sys.argv[3])  # if = 0, particles distributed in 3 dimensions, if 1, all along one straight line
 L       = N**(1./3.)
 
 # these are not used but necessary for I/O
@@ -120,44 +119,23 @@ 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))
-
-if (three_d == 0):
-    r[:,0] = radius * stheta * numpy.cos(phi)
-    r[:,1] = radius * stheta * numpy.sin(phi)
-    r[:,2] = radius * ctheta
-else:
-    r[:,0] = radius
+#r[:,0] = radius * stheta * numpy.cos(phi)
+#r[:,1] = radius * stheta * numpy.sin(phi)
+#r[:,2] = radius * ctheta
+r[:,0] = radius
 #
-
 speed  = vrot
 v      = numpy.zeros((numPart, 3))
+omega  = speed / radius
+period = 2.*math.pi/omega
+print 'period = minimum = ',min(period), ' maximum = ',max(period)
+
+omegav = omega
+if (icirc != 0):
+    omegav = omega * numpy.random.rand(N)
 
-if(three_d == 0):
-    #make N random angular momentum vectors
-    omega = numpy.random.rand(N,3)
-    #cross with r to get velocities perpendicular to displacement vector
-    v = numpy.cross(omega,r)
-    #rescale v
-    mag_v = numpy.zeros(N)
-    for i in range(N):
-        mag_v[i] = numpy.sqrt(v[i,0]**2 + v[i,1]**2 + v[i,2]**2)
-        v[i,:] /= mag_v[i]
-    v *= speed
-    if (icirc != 0):
-        for i in range(N):
-            v[i,:] *= numpy.random.rand()
-    
-else:
-    omega  = speed / radius
-    period = 2.*math.pi/omega
-    print 'period = minimum = ',min(period), ' maximum = ',max(period)
-
-    omegav = omega
-    if (icirc != 0):
-        omegav = omega * numpy.random.rand(N)
-        
-    v[:,0] = -omegav * r[:,1]
-    v[:,1] =  omegav * r[:,0]
+v[:,0] = -omegav * r[:,1]
+v[:,1] =  omegav * r[:,0]
 
 ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
 ds[()] = v

examples/IsothermalPotential/GravityOnly/plots/energy_plot.py

-import numpy as np
-import matplotlib.pyplot as plt
-import h5py as h5
-
-n_snaps = 102
-directory = '../'
-snap_name_base = "Isothermal_"
-plot_dir = "random_orbits_1d"
-plot_filename_base = "snap"
-
-v_rot = 200.
-r_0 = 100.
-
-e_kin_array = []
-e_pot_array = []
-e_total_array = []
-snap_numbers = np.linspace(0,n_snaps)
-for i in range(n_snaps):
-    snap_number = "%03d" %i
-    filename = directory + snap_name_base + snap_number + ".hdf5"
-    f = h5.File(filename)
-
-    header = f["Header"]
-    n_parts = header.attrs["NumPart_ThisFile"][1]
-    pos = np.array(f["PartType1/Coordinates"])
-    x = pos[:,0]
-    y = pos[:,1]
-    z = pos[:,2]
-    r = np.sqrt(x**2 + y**2 + z**2)
-    e_pot_particles = v_rot**2 * np.log(r/r_0)
-    e_pot_total = np.sum(e_pot_particles)
-    e_pot_array = np.append(e_pot_array,e_pot_total)
-    v = np.array(f["PartType1/Velocities"])
-    e_kin_particles = 0.5*v**2
-    e_kin_total = np.sum(e_kin_particles)
-    e_kin_array = np.append(e_kin_array,e_kin_total)
-    e_total_array = np.append(e_total_array,e_kin_total+e_pot_total)
-
-print snap_numbers.shape
-print e_kin_array.shape
-plt.plot(e_kin_array,label = "Kinetic Energy")
-plt.plot(e_pot_array,label = "Potential Energy")
-plt.plot(e_total_array,label = "Total Energy")
-plt.legend(loc = "upper right")    
-plot_filename = "%s/energy_plot.png" %plot_dir
-plt.savefig(plot_filename,format = "png")
-plt.close()

examples/IsothermalPotential/GravityOnly/plots/plot_particle_positions.py

-import h5py as h5
-import numpy as np
-import matplotlib.pyplot as plt
-
-n_snaps = 102
-directory = '../'
-snap_name_base = "Isothermal_"
-plot_directory = "random_orbits_3d"
-plot_filename_base = "snap"
-for i in range(n_snaps):
-    snap_number = "%03d" %i
-    filename = directory + snap_name_base + snap_number + ".hdf5"
-    f = h5.File(filename)
-    #find the box size
-
-    header = f["Header"]
-    box_size = header.attrs["BoxSize"][0]
-    pos = np.array(f["PartType1/Coordinates"])
-    x = pos[:,0]
-    y = pos[:,1]
-    z = pos[:,2]
-
-    #x_projection
-    plt.plot(y,z,'bo')
-    plt.xlabel("y")
-    plt.ylabel("z")
-    plt.xlim((0.0,box_size))
-    plt.ylim((0.0,box_size))
-    plot_dir = "%s/projection_x" %plot_directory
-    plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
-
-    #y_projection
-    plt.plot(x,z,'bo')
-    plt.xlabel("x")
-    plt.ylabel("z")
-    plt.xlim((0.0,box_size))
-    plt.ylim((0.0,box_size))
-    plot_dir = "%s/projection_y" %plot_directory
-    plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
-
-    #x_projection
-    plt.plot(x,y,'bo')
-    plt.xlabel("x")
-    plt.ylabel("y")
-    plt.xlim((0.0,box_size))
-    plt.ylim((0.0,box_size))
-    plot_dir = "%s/projection_z" %plot_directory
-    plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
-    
-

examples/IsothermalPotential/GravityOnly/run.sh

@@ -4,7 +4,7 @@
 if [ ! -e Isothermal.hdf5 ]
 then
     echo "Generating initial conditions for the isothermal potential box example..."
-    python makeIC.py 10000 1 1
+    python makeIC.py 1000 1
 fi
 
 ../../swift -g -t 2 isothermal.yml 2>&1 | tee output.log