Commit 3ad59a18 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch '2d_hydrodynamics' into 'master'

Hydrodynamics in 2D + New test cases

This adds the ability to run SWIFT in 1D and 2D. There is a large number of hydrodynamics tests that only exist in 2D so we will now be able to run them. This implements #145.

I have also renamed a few of the test-cases to specify the dimensionality of the problem. 

The Gresho-Chan vortex test-case can now be run and give the expected answer (provided the code is compiled with 2d hydro in const.h). The run.sh script runs the code, plots the solution and makes a movie out of everything.

See merge request !216
parents 0856bdbd 87b94377
......@@ -29,6 +29,7 @@ examples/used_parameters.yml
examples/energy.txt
examples/*/*.xmf
examples/*/*.hdf5
examples/*/*.png
examples/*/*.txt
examples/*/used_parameters.yml
examples/*/*/*.xmf
......
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_128.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
# 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).
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units).
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
......@@ -21,93 +21,83 @@
import h5py
import random
from numpy import *
import sys
# Generates a swift IC file for the Gresho Vortex in a periodic box
# Generates a swift IC file for the Gresho-Chan vortex in a periodic box
# Parameters
periodic= 1 # 1 For periodic box
factor = 3
boxSize = [ 1.0 , 1.0, 1.0/factor ]
L = 120 # Number of particles along one axis
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
rho = 1 # Gas density
rho0 = 1 # Gas density
P0 = 0. # Constant additional pressure (should have no impact on the dynamics)
fileName = "greshoVortex.hdf5"
vol = boxSize[0] * boxSize[1] * boxSize[2]
fileOutputName = "greshoVortex.hdf5"
fileGlass = "glassPlane_128.hdf5"
#---------------------------------------------------
numPart = L**3 / factor
mass = boxSize[0]*boxSize[1]*boxSize[2] * rho / numPart
#Generate particles
coords = zeros((numPart, 3))
v = zeros((numPart, 3))
m = zeros((numPart, 1))
h = zeros((numPart, 1))
u = zeros((numPart, 1))
ids = zeros((numPart, 1), dtype='L')
partId=0
for i in range(L):
for j in range(L):
for k in range(L/factor):
index = i*L*L/factor + j*L/factor + k
x = i * boxSize[0] / L + boxSize[0] / (2*L)
y = j * boxSize[0] / L + boxSize[0] / (2*L)
z = k * boxSize[0] / L + boxSize[0] / (2*L)
r2 = (x - boxSize[0] / 2)**2 + (y - boxSize[1] / 2)**2
r = sqrt(r2)
coords[index,0] = x
coords[index,1] = y
coords[index,2] = z
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[index,0] = -v_phi * (y - boxSize[0] / 2) / r
v[index,1] = v_phi * (x - boxSize[0] / 2) / r
v[index,2] = 0.
m[index] = mass
h[index] = eta * boxSize[0] / L
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[index] = P / ((gamma - 1.)*rho)
ids[index] = partId + 1
partId = partId + 1
# 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**2 / numPart
u = zeros(numPart)
v = zeros((numPart, 3))
for i in range(numPart):
x = coords[i,0]
y = coords[i,1]
z = coords[i,2]
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
file = h5py.File(fileName, 'w')
fileOutput = h5py.File(fileOutputName, 'w')
# Header
grp = file.create_group("/Header")
grp = fileOutput.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["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]
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
grp = fileOutput.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = 1
#Units
grp = file.create_group("/Units")
grp = fileOutput.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.
......@@ -115,20 +105,20 @@ 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 = fileOutput.create_group("/PartType0")
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
ds[()] = coords
ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
ds = grp.create_dataset('Masses', (numPart,1), 'f')
ds[()] = m
ds = grp.create_dataset('Masses', (numPart, 1), 'f')
ds[()] = m.reshape((numPart,1))
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h
ds[()] = h.reshape((numPart,1))
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
ds[()] = u.reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
ds[()] = ids[:]
ds[()] = ids.reshape((numPart,1))
file.close()
fileOutput.close()
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@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/>.
#
##############################################################################
# Computes the analytical solution of the Gresho-Chan vortex and plots the SPH answer
# Parameters
gas_gamma = 5./3. # Gas adiabatic index
rho0 = 1 # Gas density
P0 = 0. # Constant additional pressure (should have no impact on the dynamics)
# ---------------------------------------------------------------
# Don't touch anything after this.
# ---------------------------------------------------------------
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' : (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])
# Generate the analytic solution at this time
N = 200
R_max = 0.8
solution_r = arange(0, R_max, R_max / N)
solution_P = zeros(N)
solution_v_phi = zeros(N)
solution_v_r = zeros(N)
for i in range(N):
if solution_r[i] < 0.2:
solution_P[i] = P0 + 5. + 12.5*solution_r[i]**2
solution_v_phi[i] = 5.*solution_r[i]
elif solution_r[i] < 0.4:
solution_P[i] = P0 + 9. + 12.5*solution_r[i]**2 - 20.*solution_r[i] + 4.*log(solution_r[i]/0.2)
solution_v_phi[i] = 2. -5.*solution_r[i]
else:
solution_P[i] = P0 + 3. + 4.*log(2.)
solution_v_phi[i] = 0.
solution_rho = ones(N) * rho0
solution_s = solution_P / solution_rho**gas_gamma
solution_u = solution_P /((gas_gamma - 1.)*solution_rho)
# Read the simulation data
sim = h5py.File("gresho_%03d.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"]
pos = sim["/PartType0/Coordinates"][:,:]
x = pos[:,0] - boxSize / 2
y = pos[:, 1] - boxSize / 2
vel = sim["/PartType0/Velocities"][:,:]
r = sqrt(x**2 + y**2)
v_r = (x * vel[:,0] + y * vel[:,1]) / r
v_phi = (-y * vel[:,0] + x * vel[:,1]) / r
v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
rho = sim["/PartType0/Density"][:]
u = sim["/PartType0/InternalEnergy"][:]
S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
# Plot the interesting quantities
figure()
# Azimuthal velocity profile -----------------------------
subplot(231)
plot(r, v_phi, '.', color='r', ms=0.5)
plot(solution_r, solution_v_phi, '--', color='k', alpha=0.8, lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Azimuthal~velocity}}~v_\\phi$", labelpad=0)
xlim(0,R_max)
ylim(-0.1, 1.2)
# Radial density profile --------------------------------
subplot(232)
plot(r, rho, '.', color='r', ms=0.5)
plot(solution_r, solution_rho, '--', color='k', alpha=0.8, lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
xlim(0,R_max)
ylim(rho0-0.3, rho0 + 0.3)
#yticks([-0.2, -0.1, 0., 0.1, 0.2])
# Radial pressure profile --------------------------------
subplot(233)
plot(r, P, '.', color='r', ms=0.5)
plot(solution_r, solution_P, '--', color='k', alpha=0.8, lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Pressure}}~P$", labelpad=0)
xlim(0, R_max)
ylim(4.9 + P0, P0 + 6.1)
# Internal energy profile --------------------------------
subplot(234)
plot(r, u, '.', color='r', ms=0.5)
plot(solution_r, solution_u, '--', color='k', alpha=0.8, lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
xlim(0,R_max)
ylim(7.3, 9.1)
# Radial entropy profile --------------------------------
subplot(235)
plot(r, S, '.', color='r', ms=0.5)
plot(solution_r, solution_s, '--', color='k', alpha=0.8, lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Entropy}}~S$", labelpad=0)
xlim(0, R_max)
ylim(4.9 + P0, P0 + 6.1)
# Image --------------------------------------------------
#subplot(234)
#scatter(pos[:,0], pos[:,1], c=v_norm, cmap="PuBu", edgecolors='face', s=4, vmin=0, vmax=1)
#text(0.95, 0.95, "$|v|$", ha="right", va="top")
#xlim(0,1)
#ylim(0,1)
#xlabel("$x$", labelpad=0)
#ylabel("$y$", labelpad=0)
# Information -------------------------------------
subplot(236, frameon=False)
text(-0.49, 0.9, "Gresho-Chan vortex with $\\gamma=%.3f$ at $t=%.2f$"%(gas_gamma,time), fontsize=10)
text(-0.49, 0.8, "Background $\\rho_0=%.3f$"%rho0, fontsize=10)
text(-0.49, 0.7, "Background $P_0=%.3f$"%P0, 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([])
savefig("GreshoVortex.png", dpi=200)
#!/bin/bash
# Generate the initial conditions if they are not present.
if [ ! -e glass_128.hdf5 ]
then
echo "Fetching initial glass file for the Gresho-Chan vortex example..."
./getGlass.sh
fi
if [ ! -e greshoVortex.hdf5 ]
then
echo "Generating initial conditions for the Gresho-Chan vortex example..."
python makeIC.py
fi
# Run SWIFT
../swift -s -t 1 gresho.yml
# Plot the solution
python plotSolution.py 11
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@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 random
from numpy import *
# Computes the analytical solution of the Gresho-Chan vortex
# The script works for a given initial box and background pressure and computes the solution for any time t (The solution is constant over time).
# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
# entropic function on N radial points between r=0 and r=R_max.
# Parameters
rho0 = 1. # Background Density
P0 = 0. # Background Pressure
gamma = 5./3. # Gas polytropic index
N = 1000 # Number of radial points
R_max = 1. # Maximal radius
# ---------------------------------------------------------------
# Don't touch anything after this.
# ---------------------------------------------------------------
r = arange(0, R_max, R_max / N)
rho = ones(N)
P = zeros(N)
v = zeros(N)
u = zeros(N)
s = zeros(N)
for i in range(N):
if r[i] < 0.2:
P[i] = P0 + 5. + 12.5*r[i]**2
v[i] = 5.*r[i]
elif r[i] < 0.4:
P[i] = P0 + 9. + 12.5*r[i]**2 - 20.*r[i] + 4.*log(r[i]/0.2)
v[i] = 2. -5.*r[i]
else:
P[i] = P0 + 3. + 4.*log(2.)
v[i] = 0.
rho[i] = rho0
s[i] = P[i] / rho[i]**gamma
u[i] = P[i] /((gamma - 1.)*rho[i])
savetxt("rho.dat", column_stack((r, rho)))
savetxt("P.dat", column_stack((r, P)))
savetxt("v.dat", column_stack((r, v)))
savetxt("u.dat", column_stack((r, u)))
savetxt("s.dat", column_stack((r, s)))
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@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
import random
from numpy import *
# Generates a swift IC file containing a perturbed 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. # Pressure
gamma = 5./3. # Gas adiabatic index
pert = 0.1 # Perturbation scale (in units of the interparticle separation)
fileName = "perturbedPlane.hdf5"
#---------------------------------------------------
numPart = L**2
mass = boxSize**2 * rho / numPart
internalEnergy = P / ((gamma - 1.)*rho)
#Generate particles
coords = zeros((numPart, 3))
v = zeros((numPart, 3))
m = zeros((numPart, 1))
h = zeros((numPart, 1))
u = zeros((numPart, 1))
ids = zeros((numPart, 1), dtype='L')
for i in range(L):
for j in range(L):
index = i*L + j
x = i * boxSize / L + boxSize / (2*L) + random.random() * pert * boxSize/(2.*L)
y = j * boxSize / L + boxSize / (2*L) + random.random() * pert * boxSize/(2.*L)
z = 0
coords[index,0] = x
coords[index,1] = y
coords[index,2] = z
v[index,0] = 0.
v[index,1] = 0.
v[index,2] = 0.
m[index] = mass
h[index] = 1.23485 * boxSize / L
u[index] = internalEnergy
ids[index] = index
#--------------------------------------------------
#File
file = h5py.File(fileName, 'w')
# Header