Commit 5dad5816 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge remote-tracking branch 'origin/master' into restart-structs

parents 1dc820a7 079f0ea3
......@@ -34,13 +34,17 @@ directory. See README for run parameters.
SWIFT has been successfully built and tested with the following compilers:
- GCC 4.8.x
- GCC 4.8.x
- Intel ICC 15.0.x
- clang 3.4.x
- clang 3.4.x
More recent versions and slightly older ones should also be able to
build the software.
It has also been built with Intel and GNU C++ compilers, but that currently
requires the --disable-vec and, for Intel, --disable-compiler-warnings
configure options.
By default an attempt to choose suitable set of optimizing compiler flags
will be made, targeted for the host machine of the build. If this doesn't
work or the binaries will for another architecture then you can stop the
......@@ -61,7 +65,7 @@ You could also add some additional flags:
./configure --enable-debug --disable-optimization CFLAGS="-O2"
for instance. GCC address sanitizer flags can be included using the
for instance. GCC address sanitizer flags can be included using the
./configure --enable-sanitizer
......@@ -113,7 +117,7 @@ before you can build it.
none-standard name.
- libtool: The build system relies on libtool.
- libtool: The build system relies on libtool as well as the other autotools.
Optional Dependencies
......
......@@ -54,6 +54,9 @@ AC_USE_SYSTEM_EXTENSIONS
AX_COMPILER_VENDOR
AX_COMPILER_VERSION
# Restrict support.
AC_C_RESTRICT
# Interprocedural optimization support. Needs special handling for linking and
# archiving as well as compilation with Intels, needs to be done before
# libtool is configured (to use correct LD).
......@@ -454,6 +457,14 @@ AC_SUBST([METIS_LIBS])
AC_SUBST([METIS_INCS])
AM_CONDITIONAL([HAVEMETIS],[test -n "$METIS_LIBS"])
# METIS fixed width integer printing can require this, so define. Only needed
# for some non C99 compilers, i.e. C++ pre C++11.
AH_VERBATIM([__STDC_FORMAT_MACROS],
[/* Needed to get PRIxxx macros from stdint.h when not using C99 */
#ifndef __STDC_FORMAT_MACROS
#define __STDC_FORMAT_MACROS 1
#endif])
# Check for grackle.
have_grackle="no"
AC_ARG_WITH([grackle],
......
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9885e39 # 10^6 solar masses
UnitMass_in_cgs: 1.98848e39 # 10^6 solar masses
UnitLength_in_cgs: 3.0856776e21 # Kiloparsecs
UnitVelocity_in_cgs: 1e5 # Kilometres per second
UnitCurrent_in_cgs: 1 # Amperes
......
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9885e39 # 10^6 solar masses
UnitMass_in_cgs: 1.98848e39 # 10^6 solar masses
UnitLength_in_cgs: 3.0856776e21 # Kiloparsecs
UnitVelocity_in_cgs: 1e5 # Kilometres per second
UnitCurrent_in_cgs: 1 # Amperes
......
# 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
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 0.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-3 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: evrard # 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)
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
epsilon: 0.001 # Softening length (in internal units).
theta: 0.9
a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
r_cut: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
# Parameters related to the initial conditions
InitialConditions:
file_name: ./evrard.hdf5 # The file to read
PhysicalConstants:
G: 1.
#! /bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/evrardCollapse3D_exact.txt
################################################################################
# This file is part of SWIFT.
# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
#
# 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
from numpy import *
# Generates a swift IC file for the Evrard collapse
# Parameters
gamma = 5. / 3. # Gas adiabatic index
M = 1. # total mass of the sphere
R = 1. # radius of the sphere
u0 = 0.05 / M # initial thermal energy
fileName = "evrard.hdf5"
numPart = 100000
r = R * sqrt(random.random(numPart))
phi = 2. * pi * random.random(numPart)
cos_theta = 2. * random.random(numPart) - 1.
sin_theta = sqrt(1. - cos_theta**2)
cos_phi = cos(phi)
sin_phi = sin(phi)
pos = zeros((numPart, 3))
pos[:,0] = r * sin_theta * cos_phi
pos[:,1] = r * sin_theta * sin_phi
pos[:,2] = r * cos_theta
# shift particles to put the sphere in the centre of the box
pos += array([50. * R, 50. * R, 50. * R])
h = ones(numPart) * 2. * R / numPart**(1. / 3.)
# Generate extra arrays
v = zeros((numPart, 3))
ids = linspace(1, numPart, numPart)
m = ones(numPart) * M / numPart
u = ones(numPart) * u0
#--------------------------------------------------
#File
file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = [100. * R, 100. * R, 100. * R]
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
grp.attrs["Dimension"] = 3
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = 0
#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")
grp.create_dataset('Coordinates', data=pos, dtype='d')
grp.create_dataset('Velocities', data=v, dtype='f')
grp.create_dataset('Masses', data=m, dtype='f')
grp.create_dataset('SmoothingLength', data=h, dtype='f')
grp.create_dataset('InternalEnergy', data=u, dtype='f')
grp.create_dataset('ParticleIDs', data=ids, dtype='L')
file.close()
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
#
# 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/>.
#
################################################################################
# Compares the swift result for the 2D spherical Sod shock with a high
# resolution 2D reference result
import matplotlib
matplotlib.use("Agg")
from pylab import *
from scipy import stats
import h5py
# Parameters
gas_gamma = 5./3. # Polytropic index
rho_L = 1. # Density left state
rho_R = 0.125 # Density right state
v_L = 0. # Velocity left state
v_R = 0. # Velocity right state
P_L = 1. # Pressure left state
P_R = 0.1 # Pressure right state
# 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' : (9.90,6.45),
'figure.subplot.left' : 0.045,
'figure.subplot.right' : 0.99,
'figure.subplot.bottom' : 0.05,
'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']})
snap = int(sys.argv[1])
# Read the simulation data
sim = h5py.File("evrard_%04d.hdf5"%snap, "r")
boxSize = sim["/Header"].attrs["BoxSize"][0]
time = sim["/Header"].attrs["Time"][0]
scheme = sim["/HydroScheme"].attrs["Scheme"]
kernel = sim["/HydroScheme"].attrs["Kernel function"]
neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
eta = sim["/HydroScheme"].attrs["Kernel eta"]
git = sim["Code"].attrs["Git Revision"]
coords = sim["/PartType0/Coordinates"]
x = sqrt((coords[:,0] - 0.5 * boxSize)**2 + (coords[:,1] - 0.5 * boxSize)**2 + \
(coords[:,2] - 0.5 * boxSize)**2)
vels = sim["/PartType0/Velocities"]
v = sqrt(vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2)
u = sim["/PartType0/InternalEnergy"][:]
S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
rho = sim["/PartType0/Density"][:]
# Bin the data
x_bin_edge = logspace(-3., log10(2.), 100)
x_bin = 0.5*(x_bin_edge[1:] + x_bin_edge[:-1])
rho_bin,_,_ = stats.binned_statistic(x, rho, statistic='mean', bins=x_bin_edge)
v_bin,_,_ = stats.binned_statistic(x, v, statistic='mean', bins=x_bin_edge)
P_bin,_,_ = stats.binned_statistic(x, P, statistic='mean', bins=x_bin_edge)
S_bin,_,_ = stats.binned_statistic(x, S, statistic='mean', bins=x_bin_edge)
u_bin,_,_ = stats.binned_statistic(x, u, statistic='mean', bins=x_bin_edge)
rho2_bin,_,_ = stats.binned_statistic(x, rho**2, statistic='mean', bins=x_bin_edge)
v2_bin,_,_ = stats.binned_statistic(x, v**2, statistic='mean', bins=x_bin_edge)
P2_bin,_,_ = stats.binned_statistic(x, P**2, statistic='mean', bins=x_bin_edge)
S2_bin,_,_ = stats.binned_statistic(x, S**2, statistic='mean', bins=x_bin_edge)
u2_bin,_,_ = stats.binned_statistic(x, u**2, statistic='mean', bins=x_bin_edge)
rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
ref = loadtxt("evrardCollapse3D_exact.txt")
# Plot the interesting quantities
figure()
# Velocity profile --------------------------------
subplot(231)
semilogx(x, -v, '.', color='r', ms=0.2)
semilogx(ref[:,0], ref[:,2], "k--", alpha=0.8, lw=1.2)
errorbar(x_bin, -v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
xlim(1.e-3, 2.)
ylim(-1.7, 0.1)
# Density profile --------------------------------
subplot(232)
loglog(x, rho, '.', color='r', ms=0.2)
loglog(ref[:,0], ref[:,1], "k--", alpha=0.8, lw=1.2)
errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
xlim(1.e-3, 2.)
ylim(1.e-2, 1.e4)
# Pressure profile --------------------------------
subplot(233)
loglog(x, P, '.', color='r', ms=0.2)
loglog(ref[:,0], ref[:,3], "k--", alpha=0.8, lw=1.2)
errorbar(x_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Pressure}}~P$", labelpad=0)
xlim(1.e-3, 2.)
ylim(1.e-4, 1.e3)
# Internal energy profile -------------------------
subplot(234)
loglog(x, u, '.', color='r', ms=0.2)
loglog(ref[:,0], ref[:,3] / ref[:,1] / (gas_gamma - 1.), "k--", alpha=0.8, lw=1.2)
errorbar(x_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
xlim(1.e-3, 2.)
ylim(1.e-2, 2.)
# Entropy profile ---------------------------------
subplot(235)
semilogx(x, S, '.', color='r', ms=0.2)
semilogx(ref[:,0], ref[:,3] / ref[:,1]**gas_gamma, "k--", alpha=0.8, lw=1.2)
errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Entropy}}~S$", labelpad=0)
xlim(1.e-3, 2.)
ylim(0., 0.25)
# Information -------------------------------------
subplot(236, frameon=False)
text(-0.49, 0.9, "Evrard collapse with $\\gamma=%.3f$ in 3D\nat $t=%.2f$"%(gas_gamma,time), fontsize=10)
plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
text(-0.49, 0.4, scheme, fontsize=10)
text(-0.49, 0.3, kernel, fontsize=10)
text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
xlim(-0.5, 0.5)
ylim(0, 1)
xticks([])
yticks([])
tight_layout()
savefig("EvrardCollapse.png", dpi=200)
#!/bin/bash
# Generate the initial conditions if they are not present.
if [ ! -e evrard.hdf5 ]
then
echo "Generating initial conditions for the Evrard collapse example..."
python makeIC.py
fi
# Run SWIFT
../swift -s -G -t 4 evrard.yml 2>&1 | tee output.log
# Get the high resolution 1D reference result if not present.
if [ ! -e evrardCollapse3D_exact.txt ]
then
echo "Fetching the reference result for the Evrard collapse example..."
./getReference.sh
fi
# Plot the solution
python plotSolution.py 8
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.0856776e21 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitMass_in_cgs: 1.98848e33 # M_sun
UnitLength_in_cgs: 3.0856776e21 # kpc
UnitVelocity_in_cgs: 1e5 # km/s
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
......
......@@ -28,7 +28,7 @@ import random
# physical constants in cgs
NEWTON_GRAVITY_CGS = 6.67408e-8
SOLAR_MASS_IN_CGS = 1.9885e33
SOLAR_MASS_IN_CGS = 1.98848e33
PARSEC_IN_CGS = 3.0856776e18
# choice of units
......@@ -86,7 +86,7 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 3.0856776e21
grp.attrs["Unit mass in cgs (U_M)"] = 1.9885e33
grp.attrs["Unit mass in cgs (U_M)"] = 1.98848e33
grp.attrs["Unit time in cgs (U_t)"] = 3.0856776e16
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
......
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.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
Scheduler:
max_top_level_cells: 15
# 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).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: gresho # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1e-1 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
InitialConditions:
file_name: ./greshoVortex.hdf5 # The file to read
################################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
#
# 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
from numpy import *
# Generates a swift IC file for the Gresho-Chan vortex in a periodic box
# Parameters
gamma = 5./3. # Gas adiabatic index
rho0 = 1 # Gas density
P0 = 0. # Constant additional pressure (should have no impact on the dynamics)
fileOutputName = "greshoVortex.hdf5"
fileGlass = "glassCube_64.hdf5"
#---------------------------------------------------
# Get position and smoothing lengths from the glass
fileInput = h5py.File(fileGlass, 'r')
coords = fileInput["/PartType0/Coordinates"][:,:]
h = fileInput["/PartType0/SmoothingLength"][:]
ids = fileInput["/PartType0/ParticleIDs"][:]
boxSize = fileInput["/Header"].attrs["BoxSize"][0]
numPart = size(h)
fileInput.close()
# Now generate the rest
m = ones(numPart) * rho0 * boxSize**3 / numPart
u = zeros(numPart)
v = zeros((numPart, 3))
for i in range(numPart):
x = coords[i,0]
y = coords[i,1]
r2 = (x - boxSize / 2)**2 + (y - boxSize / 2)**2
r = sqrt(r2)
v_phi = 0.
if r < 0.2:
v_phi = 5.*r
elif r < 0.4:
v_phi = 2. - 5.*r
else:
v_phi = 0.
v[i,0] = -v_phi * (y - boxSize / 2) / r
v[i,1] = v_phi * (x - boxSize / 2) / r
v[i,2] = 0.
P = P0
if r < 0.2:
P = P + 5. + 12.5*r2
elif r < 0.4:
P = P + 9. + 12.5*r2 - 20.*r + 4.*log(r/0.2)
else:
P = P + 3. + 4.*log(2.)
u[i] = P / ((gamma - 1.)*rho0)
#File
fileOutput = h5py.File(fileOutputName, 'w')
# Header
grp = fileOutput.create_group("/Header")
grp.attrs["BoxSize"] = [boxSize, 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["NumFileOutputsPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
grp.attrs["Dimension"] = 3
#Runtime parameters
grp = fileOutput.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = 1
#Units
grp = fileOutput.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 1.
grp.attrs["Unit mass in cgs (U_M)"] = 1.