Commit 7173903b authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' into no_unnecessary_loop_in_pair

parents eb1dfd73 f9c1d350
......@@ -28,6 +28,7 @@ Valid options are:
-G Run with self-gravity
-n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop.
-s Run with SPH
-S Run with stars
-t {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
-v [12] Increase the level of verbosity
1: MPI-rank 0 writes
......
......@@ -762,6 +762,7 @@ WARN_LOGFILE =
INPUT = @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
INPUT += @top_srcdir@/src/hydro/Minimal
INPUT += @top_srcdir@/src/gravity/Default
INPUT += @top_srcdir@/src/stars/Default
INPUT += @top_srcdir@/src/riemann
INPUT += @top_srcdir@/src/potential/point_mass
INPUT += @top_srcdir@/src/cooling/const_du
......
ICs extracted from the EAGLE suite of simulations.
WARNING: The ICs are 217GB in size. They contain ~3.4G DM particles,
~3.2G gas particles and ~170M star particles
The particle distribution here is the snapshot 27 (z=0.1) of the 100Mpc
Ref-model. h- and a- factors from the original Gadget code have been
corrected for. Variables not used in a pure hydro & gravity code have
been removed.
Everything is ready to be run without cosmological integration.
The particle load of the main EAGLE simulation can be reproduced by
running these ICs on 4096 cores.
MD5 checksum of the ICs:
2301ea73e14207b541bbb04163c5269e EAGLE_ICs_100.hdf5
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 5e-2 # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
time_end: 1e-2 # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: Feedback # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
basename: eagle # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-3 # Time between statistics output
delta_time: 1e-2 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
......@@ -31,13 +31,5 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: ./Feedback.hdf5 # The file to read
file_name: ./EAGLE_ICs_100.hdf5 # The file to read
# Parameters for feedback
SN:
time: 0.001 # time the SN explodes (internal units)
energy: 1.0 # energy of the explosion (internal units)
x: 0.5 # x-position of explostion (internal units)
y: 0.5 # y-position of explostion (internal units)
z: 0.5 # z-position of explostion (internal units)
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/EAGLE_ICs_100.hdf5
#!/bin/bash
# Generate the initial conditions if they are not present.
if [ ! -e EAGLE_ICs_100.hdf5 ]
then
echo "Fetching initial conditions for the EAGLE 100Mpc example..."
./getIC.sh
fi
../swift -s -t 16 eagle_100.yml 2>&1 | tee output.log
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)
......@@ -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).
......@@ -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.
......
......@@ -24,10 +24,10 @@ 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.672e-8
NEWTON_GRAVITY_CGS = 6.67408e-8
SOLAR_MASS_IN_CGS = 1.9885e33
PARSEC_IN_CGS = 3.0856776e18
......@@ -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
periodic = 1 # 1 For periodic box
boxSize = 100. #
max_radius = boxSize / 4. # maximum radius of particles
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,25 +92,26 @@ 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)
stheta = numpy.sqrt(1.-ctheta**2)
phi = 2 * math.pi * numpy.random.rand(N)
radius = max_radius * (numpy.random.rand(numPart))**(1./3.)
print '---------------------'
print 'Radius: minimum = ',min(radius)
print 'Radius: maximum = ',max(radius)
radius = numpy.sort(radius)
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)
v = numpy.zeros((numPart, 3))
#generate particle velocities
speed = numpy.sqrt(const_G * Mass / radius)
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 = numpy.zeros((numPart, 3))
v[:,0] = -omega * r[:,1]
v[:,1] = omega * r[:,0]
......@@ -129,17 +124,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
......
#!/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
fi
rm -rf pointMass_*.hdf5
../swift -g -t 1 externalPointMass.yml 2>&1 | tee output.log
python energy_plot.py
;
; 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
base = 'Feedback'
inf = 'Feedback_005.hdf5'
blast = [5.650488e-01, 5.004371e-01, 5.010494e-01] ; location of blast
pos = h5rd(inf,'PartType0/Coordinates')
vel = h5rd(inf,'PartType0/Velocities')
rho = h5rd(inf,'PartType0/Density')
utherm = h5rd(inf,'PartType0/InternalEnergy')
; shift to centre
for ic=0,2 do pos[ic,*] = pos[ic,*] - blast[ic]
;; distance from centre
dist = fltarr(n_elements(rho))
for ic=0,2 do dist = dist + pos[ic,*]^2
dist = sqrt(dist)
; radial velocity
vr = fltarr(n_elements(rho))
for ic=0,2 do vr = vr + pos[ic,*]*vel[ic,*]
vr = vr / dist
;
end
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# 2016 Tom Theuns (tom.theuns@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
import h5py
import sys
from numpy import *
# Generates a swift IC file containing a cartesian distribution of particles
# at a constant density and pressure in a cubic box
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 1.
L = int(sys.argv[1]) # Number of particles along one axis
rho = 1. # Density
P = 1.e-6 # Pressure
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "Feedback.hdf5"
#---------------------------------------------------
numPart = L**3
mass = boxSize**3 * rho / numPart
internalEnergy = P / ((gamma - 1.)*rho)
#--------------------------------------------------
#File
file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 1.
grp.attrs["Unit mass in cgs (U_M)"] = 1.
grp.attrs["Unit time in cgs (U_t)"] = 1.
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
#Particle group
grp = file.create_group("/PartType0")
v = zeros((numPart, 3))
ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
v = zeros(1)
m = full((numPart, 1), mass)
ds = grp.create_dataset('Masses', (numPart,1), 'f')
ds[()] = m
m = zeros(1)
h = full((numPart, 1), eta * boxSize / L)
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h
h = zeros(1)
u = full((numPart, 1), internalEnergy)
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
u = zeros(1)
ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
ds[()] = ids + 1
x = ids % L;
y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2;
coords = zeros((numPart, 3))
coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
ds[()] = coords
file.close()
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,