Commit 352f3fa0 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Updated the 3D Sod shock example to use the new glass files. Added script to...

Updated the 3D Sod shock example to use the new glass files. Added script to plot solution. Removed old glass files.
parent 839a3bcb
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
# 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 2D Sedov blast wave.
# The script works for a given initial box and dumped energy and computes the solution at a later time t.
# Parameters
rho_0 = 1. # Background Density
P_0 = 1.e-6 # Background Pressure
E_0 = 1. # Energy of the explosion
gas_gamma = 5./3. # Gas polytropic index
# ---------------------------------------------------------------
# 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])
# Read the simulation data
sim = h5py.File("sedov_%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
z = pos[:,2] - boxSize / 2
vel = sim["/PartType0/Velocities"][:,:]
r = sqrt(x**2 + y**2 + z**2)
v_r = (x * vel[:,0] + y * vel[:,1] + z * vel[:,2]) / r
u = sim["/PartType0/InternalEnergy"][:]
S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
rho = sim["/PartType0/Density"][:]
# Now, work our the solution....
from scipy.special import gamma as Gamma
from numpy import *
def calc_a(g,nu=3):
"""
exponents of the polynomials of the sedov solution
g - the polytropic gamma
nu - the dimension
"""
a = [0]*8
a[0] = 2.0 / (nu + 2)
a[2] = (1-g) / (2*(g-1) + nu)
a[3] = nu / (2*(g-1) + nu)
a[5] = 2 / (g-2)
a[6] = g / (2*(g-1) + nu)
a[1] = (((nu+2)*g)/(2.0+nu*(g-1.0)) ) * ( (2.0*nu*(2.0-g))/(g*(nu+2.0)**2) - a[2])
a[4] = a[1]*(nu+2) / (2-g)
a[7] = (2 + nu*(g-1))*a[1]/(nu*(2-g))
return a
def calc_beta(v, g, nu=3):
"""
beta values for the sedov solution (coefficients of the polynomials of the similarity variables)
v - the similarity variable
g - the polytropic gamma
nu- the dimension
"""
beta = (nu+2) * (g+1) * array((0.25, (g/(g-1))*0.5,
-(2 + nu*(g-1))/2.0 / ((nu+2)*(g+1) -2*(2 + nu*(g-1))),
-0.5/(g-1)), dtype=float64)
beta = outer(beta, v)
beta += (g+1) * array((0.0, -1.0/(g-1),
(nu+2) / ((nu+2)*(g+1) -2.0*(2 + nu*(g-1))),
1.0/(g-1)), dtype=float64).reshape((4,1))
return beta
def sedov(t, E0, rho0, g, n=1000, nu=3):
"""
solve the sedov problem
t - the time
E0 - the initial energy
rho0 - the initial density
n - number of points (10000)
nu - the dimension
g - the polytropic gas gamma
"""
# the similarity variable
v_min = 2.0 / ((nu + 2) * g)
v_max = 4.0 / ((nu + 2) * (g + 1))
v = v_min + arange(n) * (v_max - v_min) / (n - 1.0)
a = calc_a(g, nu)
beta = calc_beta(v, g=g, nu=nu)
lbeta = log(beta)
r = exp(-a[0] * lbeta[0] - a[2] * lbeta[1] - a[1] * lbeta[2])
rho = ((g + 1.0) / (g - 1.0)) * exp(a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2])
p = exp(nu * a[0] * lbeta[0] + (a[5] + 1) * lbeta[3] + (a[4] - 2 * a[1]) * lbeta[2])
u = beta[0] * r * 4.0 / ((g + 1) * (nu + 2))
p *= 8.0 / ((g + 1) * (nu + 2) * (nu + 2))
# we have to take extra care at v=v_min, since this can be a special point.
# It is not a singularity, however, the gradients of our variables (wrt v) are.
# r -> 0, u -> 0, rho -> 0, p-> constant
u[0] = 0.0; rho[0] = 0.0; r[0] = 0.0; p[0] = p[1]
# volume of an n-sphere
vol = (pi ** (nu / 2.0) / Gamma(nu / 2.0 + 1)) * power(r, nu)
# note we choose to evaluate the integral in this way because the
# volumes of the first few elements (i.e near v=vmin) are shrinking
# very slowly, so we dramatically improve the error convergence by
# finding the volumes exactly. This is most important for the
# pressure integral, as this is on the order of the volume.
# (dimensionless) energy of the model solution
de = rho * u * u * 0.5 + p / (g - 1)
# integrate (trapezium rule)
q = inner(de[1:] + de[:-1], diff(vol)) * 0.5
# the factor to convert to this particular problem
fac = (q * (t ** nu) * rho0 / E0) ** (-1.0 / (nu + 2))
# shock speed
shock_speed = fac * (2.0 / (nu + 2))
rho_s = ((g + 1) / (g - 1)) * rho0
r_s = shock_speed * t * (nu + 2) / 2.0
p_s = (2.0 * rho0 * shock_speed * shock_speed) / (g + 1)
u_s = (2.0 * shock_speed) / (g + 1)
r *= fac * t
u *= fac
p *= fac * fac * rho0
rho *= rho0
return r, p, rho, u, r_s, p_s, rho_s, u_s, shock_speed
# The main properties of the solution
r_s, P_s, rho_s, v_s, r_shock, _, _, _, _ = sedov(time, E_0, rho_0, gas_gamma, 1000, 3)
# Append points for after the shock
r_s = np.insert(r_s, np.size(r_s), [r_shock, r_shock*1.5])
rho_s = np.insert(rho_s, np.size(rho_s), [rho_0, rho_0])
P_s = np.insert(P_s, np.size(P_s), [P_0, P_0])
v_s = np.insert(v_s, np.size(v_s), [0, 0])
# Additional arrays
u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy
s_s = P_s / rho_s**gas_gamma # entropic function
# Plot the interesting quantities
figure()
# Velocity profile --------------------------------
subplot(231)
plot(r, v_r, '.', color='r', ms=1.)
plot(r_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Radial~velocity}}~v_r$", labelpad=0)
xlim(0, 1.3 * r_shock)
ylim(-0.2, 3.8)
# Density profile --------------------------------
subplot(232)
plot(r, rho, '.', color='r', ms=1.)
plot(r_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Density}}~\\rho$", labelpad=2)
xlim(0, 1.3 * r_shock)
ylim(-0.2, 5.2)
# Pressure profile --------------------------------
subplot(233)
plot(r, P, '.', color='r', ms=1.)
plot(r_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Pressure}}~P$", labelpad=0)
xlim(0, 1.3 * r_shock)
ylim(-1, 12.5)
# Internal energy profile -------------------------
subplot(234)
plot(r, u, '.', color='r', ms=1.)
plot(r_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
xlim(0, 1.3 * r_shock)
ylim(-2, 22)
# Entropy profile ---------------------------------
subplot(235)
plot(r, S, '.', color='r', ms=1.)
plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Entropy}}~S$", labelpad=0)
xlim(0, 1.3 * r_shock)
ylim(-5, 50)
# Information -------------------------------------
subplot(236, frameon=False)
text(-0.49, 0.9, "Sedov blast with $\\gamma=%.3f$ in 2D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
text(-0.49, 0.8, "Background $\\rho_0=%.2f$"%(rho_0), fontsize=10)
text(-0.49, 0.7, "Energy injected $E_0=%.2f$"%(E_0), 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("Sedov.png", dpi=200)
......@@ -22,8 +22,6 @@
# Generates the analytical solution for the Sod shock test case
# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
# 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 points between x_min and x_max.
# This follows the solution given in (Toro, 2009)
......
......@@ -22,8 +22,6 @@
# Generates the analytical solution for the Sod shock test case
# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
# 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 points between x_min and x_max.
# This follows the solution given in (Toro, 2009)
......
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# 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
......@@ -21,84 +20,75 @@
import h5py
from numpy import *
# Generates a swift IC file for the Sod Shock in a periodic box
# Generates a swift IC file for the 3D Sod Shock in a periodic box
# Parameters
periodic= 1 # 1 For periodic box
factor = 8
boxSize = [ 1.0 , 1.0/factor , 1.0/factor ]
L = 100 # Number of particles along one axis
P1 = 1. # Pressure left state
P2 = 0.1795 # Pressure right state
gamma = 5./3. # Gas adiabatic index
gamma = 5./3. # Gas adiabatic index
x_min = -1.
x_max = 1.
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
fileName = "sodShock.hdf5"
vol = boxSize[0] * boxSize[1] * boxSize[2]
#---------------------------------------------------
#Read in high density glass
# glass1 = h5py.File("../Glass/glass_200000.hdf5")
glass1 = h5py.File("glass_001.hdf5")
pos1 = glass1["/PartType0/Coordinates"][:,:]
pos1 = pos1 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
glass_h1 = glass1["/PartType0/SmoothingLength"][:] / factor
#Read in high density glass
# glass2 = h5py.File("../Glass/glass_50000.hdf5")
glass2 = h5py.File("glass_002.hdf5")
pos2 = glass2["/PartType0/Coordinates"][:,:]
pos2 = pos2 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
glass_h2 = glass2["/PartType0/SmoothingLength"][:] / factor
#Generate high density region
rho1 = 1.
coord1 = append(pos1, pos1 + [0.125, 0, 0], 0)
coord1 = append(coord1, coord1 + [0.25, 0, 0], 0)
# coord1 = append(pos1, pos1 + [0, 0.5, 0], 0)
# coord1 = append(coord1, pos1 + [0, 0, 0.5], 0)
# coord1 = append(coord1, pos1 + [0, 0.5, 0.5], 0)
N1 = size(coord1)/3
v1 = zeros((N1, 3))
u1 = ones(N1) * P1 / ((gamma - 1.) * rho1)
m1 = ones(N1) * vol * 0.5 * rho1 / N1
h1 = append(glass_h1, glass_h1, 0)
h1 = append(h1, h1, 0)
#Generate low density region
rho2 = 0.25
coord2 = append(pos2, pos2 + [0.125, 0, 0], 0)
coord2 = append(coord2, coord2 + [0.25, 0, 0], 0)
# coord2 = append(pos2, pos2 + [0, 0.5, 0], 0)
# coord2 = append(coord2, pos2 + [0, 0, 0.5], 0)
# coord2 = append(coord2, pos2 + [0, 0.5, 0.5], 0)
N2 = size(coord2)/3
v2 = zeros((N2, 3))
u2 = ones(N2) * P2 / ((gamma - 1.) * rho2)
m2 = ones(N2) * vol * 0.5 * rho2 / N2
h2 = append(glass_h2, glass_h2, 0)
h2 = append(h2, h2, 0)
#Merge arrays
numPart = N1 + N2
coords = append(coord1, coord2+[0.5, 0., 0.], 0)
v = append(v1, v2,0)
h = append(h1, h2,0)
u = append(u1, u2,0)
m = append(m1, m2,0)
ids = zeros(numPart, dtype='L')
for i in range(1, numPart+1):
ids[i-1] = i
#Final operation since we come from Gadget-2 cubic spline ICs
h /= 1.825752
boxSize = (x_max - x_min)
glass_L = h5py.File("glassCube_64.hdf5", "r")
glass_R = h5py.File("glassCube_32.hdf5", "r")
pos_L = glass_L["/PartType0/Coordinates"][:,:] * 0.5
pos_R = glass_R["/PartType0/Coordinates"][:,:] * 0.5
h_L = glass_L["/PartType0/SmoothingLength"][:] * 0.5
h_R = glass_R["/PartType0/SmoothingLength"][:] * 0.5
# Merge things
aa = pos_L - array([0.5, 0., 0.])
pos_LL = append(pos_L, pos_L + array([0.5, 0., 0.]), axis=0)
pos_RR = append(pos_R, pos_R + array([0.5, 0., 0.]), axis=0)
pos = append(pos_LL - array([1.0, 0., 0.]), pos_RR, axis=0)
h_LL = append(h_L, h_L)
h_RR = append(h_R, h_R)
h = append(h_LL, h_RR)
numPart_L = size(h_LL)
numPart_R = size(h_RR)
numPart = size(h)
vol_L = 0.25
vol_R = 0.25
# Generate extra arrays
v = zeros((numPart, 3))
ids = linspace(1, numPart, numPart)
m = zeros(numPart)
u = zeros(numPart)
for i in range(numPart):
x = pos[i,0]
if x < 0: #left
u[i] = P_L / (rho_L * (gamma - 1.))
m[i] = rho_L * vol_L / numPart_L
v[i,0] = v_L
else: #right
u[i] = P_R / (rho_R * (gamma - 1.))
m[i] = rho_R * vol_R / numPart_R
v[i,0] = v_R
# Shift particles
pos[:,0] -= x_min
#File
file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["BoxSize"] = [boxSize, 0.5, 0.5]
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]
......@@ -109,7 +99,7 @@ grp.attrs["Flag_Entropy_ICs"] = 0
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
grp.attrs["PeriodicBoundariesOn"] = 1
#Units
grp = file.create_group("/Units")
......@@ -121,7 +111,7 @@ grp.attrs["Unit temperature in cgs (U_T)"] = 1.
#Particle group
grp = file.create_group("/PartType0")
grp.create_dataset('Coordinates', data=coords, dtype='d')
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')
......@@ -130,5 +120,3 @@ grp.create_dataset('ParticleIDs', data=ids, dtype='L')
file.close()
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# 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
......@@ -18,61 +17,105 @@
#
##############################################################################
import random
from numpy import *
# Computes the analytical solution of the Sod shock and plots the SPH answer
# Generates the analytical solution for the Sod shock test case
# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
# 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 points between x_min and x_max.
# This follows the solution given in (Toro, 2009)
# Parameters
rho_L = 1
P_L = 1
v_L = 0.
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
rho_R = 0.25
P_R = 0.1795
v_R = 0.
gamma = 5./3. # Polytropic index
import matplotlib
matplotlib.use("Agg")
from pylab import *
import h5py
t = 0.12 # Time of the evolution
# 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("sodShock_%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"]
x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0]
u = sim["/PartType0/InternalEnergy"][:]
S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
rho = sim["/PartType0/Density"][:]
N = 1000 # Number of points
x_min = -0.25
x_max = 0.25
x_min = -1.
x_max = 1.
x += x_min
# ---------------------------------------------------------------
# Don't touch anything after this.
# ---------------------------------------------------------------
c_L = sqrt(gamma * P_L / rho_L) # Speed of the rarefaction wave
c_R = sqrt(gamma * P_R / rho_R) # Speed of the shock front
c_L = sqrt(gas_gamma * P_L / rho_L) # Speed of the rarefaction wave
c_R = sqrt(gas_gamma * P_R / rho_R) # Speed of the shock front
# Helpful variable
Gama = (gamma - 1.) / (gamma + 1.)