Also updated external point example to python scripts.

examples/ExternalPointMass/energy_plot.py

+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (3.15,3.15),
+'figure.subplot.left'    : 0.145,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.11,
+'figure.subplot.top'     : 0.99,
+'figure.subplot.wspace'  : 0.15,
+'figure.subplot.hspace'  : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+import numpy as np
+import h5py as h5
+import sys
+
+# File containing the total energy
+stats_filename = "./energy.txt"
+
+# First snapshot
+snap_filename = "pointMass_000.hdf5"
+f = h5.File(snap_filename,'r')
+
+# Read the units parameters from the snapshot
+units = f["InternalCodeUnits"]
+unit_mass = units.attrs["Unit mass in cgs (U_M)"]
+unit_length = units.attrs["Unit length in cgs (U_L)"]
+unit_time = units.attrs["Unit time in cgs (U_t)"]
+
+G = 6.67408e-8 * unit_mass * unit_time**2 / unit_length**3
+
+# Read the header
+header = f["Header"]
+box_size = float(header.attrs["BoxSize"][0])
+
+# Read the properties of the potential
+parameters = f["Parameters"]
+mass = float(parameters.attrs["PointMassPotential:mass"])
+centre = [box_size/2, box_size/2, box_size/2]
+f.close()
+
+# Read the statistics summary
+file_energy = np.loadtxt("energy.txt")
+time_stats = file_energy[:,0]
+E_kin_stats = file_energy[:,3]
+E_pot_stats = file_energy[:,5]
+E_tot_stats = E_kin_stats + E_pot_stats
+
+# Read the snapshots
+time_snap = np.zeros(402)
+E_kin_snap = np.zeros(402)
+E_pot_snap = np.zeros(402)
+E_tot_snap = np.zeros(402)
+Lz_snap = np.zeros(402)
+
+# Read all the particles from the snapshots
+for i in range(402):
+    snap_filename = "pointMass_%0.3d.hdf5"%i
+    f = h5.File(snap_filename,'r')
+
+    pos_x = f["PartType1/Coordinates"][:,0]
+    pos_y = f["PartType1/Coordinates"][:,1]
+    pos_z = f["PartType1/Coordinates"][:,2]
+    vel_x = f["PartType1/Velocities"][:,0]
+    vel_y = f["PartType1/Velocities"][:,1]
+    vel_z = f["PartType1/Velocities"][:,2]
+    m = f["/PartType1/Masses"][:]
+    
+    r = np.sqrt((pos_x[:] - centre[0])**2 + (pos_y[:] - centre[1])**2 + (pos_z[:] - centre[2])**2)
+    Lz = (pos_x[:] - centre[0]) * vel_y[:] - (pos_y[:] - centre[1]) * vel_x[:]
+
+    time_snap[i] = f["Header"].attrs["Time"]
+    E_kin_snap[i] = np.sum(0.5 * m * (vel_x[:]**2 + vel_y[:]**2 + vel_z[:]**2))
+    E_pot_snap[i] = np.sum(-mass * m * G / r)
+    E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
+    Lz_snap[i] = np.sum(Lz)
+
+# Plot energy evolution
+figure()
+plot(time_stats, E_kin_stats, "r-", lw=0.5, label="Kinetic energy")
+plot(time_stats, E_pot_stats, "g-", lw=0.5, label="Potential energy")
+plot(time_stats, E_tot_stats, "k-", lw=0.5, label="Total energy")
+
+plot(time_snap[::10], E_kin_snap[::10], "rD", lw=0.5, ms=2)
+plot(time_snap[::10], E_pot_snap[::10], "gD", lw=0.5, ms=2)
+plot(time_snap[::10], E_tot_snap[::10], "kD", lw=0.5, ms=2)
+
+legend(loc="center right", fontsize=8, frameon=False, handlelength=3, ncol=1)
+xlabel("${\\rm{Time}}$", labelpad=0)
+ylabel("${\\rm{Energy}}$",labelpad=0)
+xlim(0, 8)
+
+savefig("energy.png", dpi=200)
+
+# Plot angular momentum evolution
+figure()
+plot(time_snap, Lz_snap, "k-", lw=0.5, ms=2)
+
+xlabel("${\\rm{Time}}$", labelpad=0)
+ylabel("${\\rm{Angular~momentum}}$",labelpad=0)
+xlim(0, 8)
+
+savefig("angular_momentum.png", dpi=200)
+
+

examples/ExternalPointMass/externalPointMass.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:   1.    # The end time of the simulation (in internal units).
+  time_end:   8.    # The end time of the simulation (in internal units).
   dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
 

examples/ExternalPointMass/makeIC.py

@@ -27,7 +27,7 @@ import random
 # Generates a random distriution of particles, for motion in an external potnetial centred at (0,0,0)
 
 # 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
 

examples/ExternalPointMass/run.sh

@@ -9,3 +9,5 @@ fi
 
 rm -rf pointMass_*.hdf5
 ../swift -g -t 1 externalPointMass.yml 2>&1 | tee output.log
+
+python energy_plot.py

examples/ExternalPointMass/test.pro

-;
-;  test energy / angular momentum conservation of test problem
-;
-@physunits
-
-indir    = '/gpfs/data/tt/Codes/Swift-git/swiftsim/examples/'
-basefile = 'output_'
-nfiles   = 657
-nfollow  = 100 ; number of particles to follow
-eout     = fltarr(nfollow, nfiles)
-ekin     = fltarr(nfollow, nfiles)
-epot     = fltarr(nfollow, nfiles)
-tout     = fltarr(nfiles)
-; set properties of potential
-uL  = 1e3 * phys.pc             ; unit of length
-uM  = phys.msun                 ; unit of mass
-uV  = 1d5                       ; unit of velocity
-
-; derived units
-constG   = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
-pcentre  = [50.,50.,50.] * 1d3 * pc / uL
-mextern  = 1d10 * msun / uM
-;
-;
-;
-ifile  = 0
-for ifile=0,nfiles-1 do begin
-;for ifile=0,3 do begin
-   inf    = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
-   time   = h5ra(inf, 'Header','Time')
-   p      = h5rd(inf,'PartType1/Coordinates')
-   v      = h5rd(inf,'PartType1/Velocities')
-   id     = h5rd(inf,'PartType1/ParticleIDs')
-   indx   = sort(id)
-;
-   id     = id[indx]
-   for ic=0,2 do begin
-      tmp = reform(p[ic,*]) & p[ic,*] = tmp[indx]
-      tmp = reform(v[ic,*]) & v[ic,*] = tmp[indx]
-   endfor
-; calculate energy
-   dd  = size(p,/dimen) & npart = dd[1]
-   ener = fltarr(npart)
-   dr   = fltarr(npart) & dv = dr
-   for ic=0,2 do dr[*] = dr[*] + (p[ic,*]-pcentre[ic])^2
-   for ic=0,2 do dv[*] = dv[*] + v[ic,*]^2
-   dr = sqrt(dr)
-;   print,'time = ',time,p[0,0],v[0,0],id[0]
-   ek   = 0.5 * dv
-   ep   = - constG * mextern / dr
-   ener = ek + ep
-   tout(ifile) = time
-   eout(*,ifile) = ener[0:nfollow-1]
-   ekin(*,ifile) = ek[0:nfollow-1]
-   epot(*,ifile) = ep[0:nfollow-1]
-endfor
-
-; calculate relative energy change
-de = 0.0 * eout
-for ifile=1, nfiles -1 do de[*,ifile] = (eout[*,ifile]-eout[*,0])/eout[*,0]
-
-
-end
-
-

examples/Makefile.am

@@ -63,7 +63,7 @@ EXTRA_DIST = BigCosmoVolume/makeIC.py \
 	     EAGLE_12/eagle_12.yml EAGLE_12/getIC.sh EAGLE_12/README EAGLE_12/run.sh \
 	     EAGLE_25/eagle_25.yml EAGLE_25/getIC.sh EAGLE_25/README EAGLE_25/run.sh \
 	     EAGLE_50/eagle_50.yml EAGLE_50/getIC.sh EAGLE_50/README EAGLE_50/run.sh \
-	     ExternalPointMass/externalPointMass.yml ExternalPointMass/makeIC.py ExternalPointMass/run.sh ExternalPointMass/test.pro \
+	     ExternalPointMass/externalPointMass.yml ExternalPointMass/makeIC.py ExternalPointMass/run.sh ExternalPointMass/energy_plot.py \
 	     GreshoVortex_2D/getGlass.sh GreshoVortex_2D/gresho.yml GreshoVortex_2D/makeIC.py GreshoVortex_2D/plotSolution.py GreshoVortex_2D/run.sh \
 	     HydrostaticHalo/README HydrostaticHalo/hydrostatic.yml HydrostaticHalo/makeIC.py HydrostaticHalo/run.sh \
 	     HydrostaticHalo/density_profile.py HydrostaticHalo/velocity_profile.py HydrostaticHalo/internal_energy_profile.py HydrostaticHalo/test_energy_conservation.py \