diff --git a/examples/ExternalPointMass/energy_plot.py b/examples/ExternalPointMass/energy_plot.py new file mode 100644 index 0000000000000000000000000000000000000000..a75fcb835d33b3695170aab822092556f12db7d1 --- /dev/null +++ b/examples/ExternalPointMass/energy_plot.py @@ -0,0 +1,121 @@ +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) + + diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml index 621a66bbc39838ac8d3d8a8a3992b2a7be3157a8..f77b4632c67115559ab2d45083eb7778d21903d8 100644 --- a/examples/ExternalPointMass/externalPointMass.yml +++ b/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). diff --git a/examples/ExternalPointMass/makeIC.py b/examples/ExternalPointMass/makeIC.py index 326183398933c88d7348e72e00343064b3e3a64c..a324220e183baed4ac5aaaf2d71cf5b6ea15c1ce 100644 --- a/examples/ExternalPointMass/makeIC.py +++ b/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 diff --git a/examples/ExternalPointMass/run.sh b/examples/ExternalPointMass/run.sh index 8b745746e40341dc019005912fe1e6a75e0e87ac..44daaaaf5b92c2923527daf256a61fbb4403b73f 100755 --- a/examples/ExternalPointMass/run.sh +++ b/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 diff --git a/examples/ExternalPointMass/test.pro b/examples/ExternalPointMass/test.pro deleted file mode 100644 index 21c10e9d27daa45b085c6a659ba3cf7260f017fb..0000000000000000000000000000000000000000 --- a/examples/ExternalPointMass/test.pro +++ /dev/null @@ -1,65 +0,0 @@ -; -; 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 - - diff --git a/examples/Makefile.am b/examples/Makefile.am index 58dcc5ac2b45e567d3023303d2a588551a5cfcfd..dd13fb7eb4b82fbbfbb1ae450e20d01b13f2a455 100644 --- a/examples/Makefile.am +++ b/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 \