Skip to content
GitLab
Commit dea2e02d authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Updated the Sedov-Taylor blast 3D example

parent d5f913b8
Hide whitespace changes
Inline Side-by-side
examples/SedovBlast_2D/makeIC.py
View file @ dea2e02d
......@@ -18,7 +18,6 @@
##############################################################################
import h5py
import random
from numpy import *
# Generates a swift IC file for the Sedov blast test in a periodic cubic box
......
examples/SedovBlast_3D/makeIC.py
View file @ dea2e02d
###############################################################################
# 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
......@@ -19,65 +18,42 @@
##############################################################################
import h5py
import random
from numpy import *
# Generates a swift IC file for the Sedov blast test in a periodic cubic box
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 10.
L = 101 # Number of particles along one axis
rho = 1. # Density
P = 1.e-5 # Pressure
E0= 1.e2 # Energy of the explosion
pert = 0.1
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
gamma = 5./3. # Gas adiabatic index
rho0 = 1. # Background density
P0 = 1.e-6 # Background pressure
E0= 1. # Energy of the explosion
N_inject = 15 # Number of particles in which to inject energy
fileName = "sedov.hdf5"
#---------------------------------------------------
numPart = L**3
mass = boxSize**3 * rho / numPart
internalEnergy = P / ((gamma - 1.)*rho)
if L%2 == 0:
print "Number of particles along each dimension must be odd."
exit()
#Generate particles
coords = zeros((numPart, 3))
v = zeros((numPart, 3))
m = zeros(numPart)
h = zeros(numPart)
u = zeros(numPart)
ids = zeros(numPart, dtype='L')
for i in range(L):
for j in range(L):
for k in range(L):
index = i*L*L + j*L + k
x = i * boxSize / L + boxSize / (2*L)
y = j * boxSize / L + boxSize / (2*L)
z = k * boxSize / L + boxSize / (2*L)
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] = eta * boxSize / L
u[index] = internalEnergy
ids[index] = index + 1
if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L:
u[index] = u[index] + E0 / (33. * mass)
print "Particle " , index , " set to detonate."
coords[index,0] = x + random.random() * pert * boxSize/(2.*L)
coords[index,1] = y + random.random() * pert * boxSize/(2.*L)
coords[index,2] = z + random.random() * pert * boxSize/(2.*L)
glass = h5py.File("glassCube_64.hdf5", "r")
# Read particle positions and h from the glass
pos = glass["/PartType0/Coordinates"][:,:]
h = glass["/PartType0/SmoothingLength"][:] * 0.3
numPart = size(h)
vol = 1.
# Generate extra arrays
v = zeros((numPart, 3))
ids = linspace(1, numPart, numPart)
m = zeros(numPart)
u = zeros(numPart)
r = zeros(numPart)
r = sqrt((pos[:,0] - 0.5)**2 + (pos[:,1] - 0.5)**2 + (pos[:,2] - 0.5)**2)
m[:] = rho0 * vol / numPart
u[:] = P0 / (rho0 * (gamma - 1))
# Make the central particles detonate
index = argsort(r)
u[index[0:N_inject]] = E0 / (N_inject * m[0])
#--------------------------------------------------
......@@ -86,7 +62,7 @@ file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["BoxSize"] = [1., 1., 1.]
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]
......@@ -97,7 +73,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")
......@@ -109,7 +85,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')
......
examples/SedovBlast_3D/makeIC_fcc.py deleted 100644 → 0
View file @ d5f913b8
###############################################################################
# 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 h5py
import random
from numpy import *
# Generates a swift IC file for the Sedov blast test in a periodic cubic box
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 5.
L = 64 # Number of particles boxes along one axis
rho = 1. # Density
P = 1.e-5 # Pressure
E0= 1.e2 # Energy of the explosion
pert = 0.025
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "sedov.hdf5"
#---------------------------------------------------
numPart = 4*(L**3)
mass = boxSize**3 * rho / numPart
internalEnergy = P / ((gamma - 1.)*rho)
off = array( [ [ 0.0 , 0.0 , 0.0 ] , [ 0.0 , 0.5 , 0.5 ] , [ 0.5 , 0.0 , 0.5 ] , [ 0.5 , 0.5 , 0.0 ] ] );
hbox = boxSize / L
# if L%2 == 0:
# print "Number of particles along each dimension must be odd."
# exit()
#Generate particles
coords = zeros((numPart, 3))
v = zeros((numPart, 3))
m = zeros(numPart)
h = zeros(numPart)
u = zeros(numPart)
ids = zeros(numPart, dtype='L')
for i in range(L):
for j in range(L):
for k in range(L):
x = (i + 0.25) * hbox
y = (j + 0.25) * hbox
z = (k + 0.25) * hbox
for ell in range(4):
index = 4*(i*L*L + j*L + k) + ell
coords[index,0] = x + off[ell,0] * hbox
coords[index,1] = y + off[ell,1] * hbox
coords[index,2] = z + off[ell,2] * hbox
v[index,0] = 0.
v[index,1] = 0.
v[index,2] = 0.
m[index] = mass
h[index] = eta * hbox
u[index] = internalEnergy
ids[index] = index + 1
if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox:
u[index] = u[index] + E0 / (28. * mass)
print "Particle " , index , " set to detonate."
coords[index,0] += random.random() * pert * hbox
coords[index,1] += random.random() * pert * hbox
coords[index,2] += random.random() * pert * hbox
#--------------------------------------------------
#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")
grp.create_dataset('Coordinates', data=coords, 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()
examples/SedovBlast_3D/profile.py deleted 100644 → 0
View file @ d5f913b8
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
#
# 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 numpy as np
import h5py
import matplotlib
matplotlib.use('Agg')
import pylab as pl
import glob
import sedov
# plot the radial density profile for all snapshots in the current directory,
# as well as the analytical solution of Sedov
for f in sorted(glob.glob("output*.hdf5")):
file = h5py.File(f, "r")
t = file["/Header"].attrs["Time"]
coords = np.array(file["/PartType0/Coordinates"])
rho = np.array(file["/PartType0/Density"])
radius = np.sqrt( (coords[:,0]-5.)**2 + (coords[:,1]-5.)**2 + \
(coords[:,2]-5.)**2 )
r_theory, rho_theory = sedov.get_analytical_solution(100., 5./3., 3, t,
max(radius))
pl.plot(r_theory, rho_theory, "r-")
pl.plot(radius, rho, "k.")
pl.savefig("{name}.png".format(name = f[:-5]))
pl.close()
examples/SedovBlast_3D/rdf.py deleted 100644 → 0
View file @ d5f913b8
###############################################################################
# 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 h5py
import random
import sys
import math
from numpy import *
# Reads the HDF5 output of SWIFT and generates a radial density profile
# of the different physical values.
# Input values?
if len(sys.argv) < 3 :
print "Usage: " , sys.argv[0] , " <filename> <nr. bins>"
exit()
# Get the input arguments
fileName = sys.argv[1];
nr_bins = int( sys.argv[2] );
# Open the file
fileName = sys.argv[1];
file = h5py.File( fileName , 'r' )
# Get the space dimensions.
grp = file[ "/Header" ]
boxsize = grp.attrs[ 'BoxSize' ]
if len( boxsize ) > 0:
boxsize = boxsize[0]
# Get the particle data
grp = file.get( '/PartType0' )
ds = grp.get( 'Coordinates' )
coords = ds[...]
ds = grp.get( 'Velocities' )
v = ds[...]
ds = grp.get( 'Masses' )
m = ds[...]
ds = grp.get( 'SmoothingLength' )
h = ds[...]
ds = grp.get( 'InternalEnergy' )
u = ds[...]
ds = grp.get( 'ParticleIDs' )
ids = ds[...]
ds = grp.get( 'Density' )
rho = ds[...]
# Set the origin to be the middle of the box
origin = [ boxsize / 2 , boxsize / 2 , boxsize / 2 ]
# Get the maximum radius
r_max = boxsize / 2
# Init the bins
nr_parts = coords.shape[0]
bins_v = zeros( nr_bins )
bins_m = zeros( nr_bins )
bins_h = zeros( nr_bins )
bins_u = zeros( nr_bins )
bins_rho = zeros( nr_bins )
bins_count = zeros( nr_bins )
bins_P = zeros( nr_bins )
# Loop over the particles and fill the bins.
for i in range( nr_parts ):
# Get the box index.
r = coords[i] - origin
r = sqrt( r[0]*r[0] + r[1]*r[1] + r[2]*r[2] )
ind = floor( r / r_max * nr_bins )
# Update the bins
if ( ind < nr_bins ):
bins_count[ind] += 1
bins_v[ind] += sqrt( v[i,0]*v[i,0] + v[i,1]*v[i,1] + v[i,2]*v[i,2] )
bins_m[ind] += m[i]
bins_h[ind] += h[i]
bins_u[ind] += u[i]
bins_rho[ind] += rho[i]
bins_P[ind] += (2.0/3)*u[i]*rho[i]
# Loop over the bins and dump them
print "# bucket left right count v m h u rho"
for i in range( nr_bins ):
# Normalize by the bin volume.
r_left = r_max * i / nr_bins
r_right = r_max * (i+1) / nr_bins
vol = 4/3*math.pi*(r_right*r_right*r_right - r_left*r_left*r_left)
ivol = 1.0 / vol
print "%i %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e" % \
( i , r_left , r_right , \
bins_count[i] * ivol , \
bins_v[i] / bins_count[i] , \
bins_m[i] * ivol , \
bins_h[i] / bins_count[i] , \
bins_u[i] / bins_count[i] , \
bins_rho[i] / bins_count[i] ,
bins_P[i] / bins_count[i] )
examples/SedovBlast_3D/run.sh
View file @ dea2e02d
#!/bin/bash
# Generate the initial conditions if they are not present.
# Generate the initial conditions if they are not present.
if [ ! -e glassCube_64.hdf5 ]
then
echo "Fetching initial glass file for the Sedov blast example..."
./getGlass.sh
fi
if [ ! -e sedov.hdf5 ]
then
echo "Generating initial conditions for the SedovBlast example..."
python makeIC_fcc.py
echo "Generating initial conditions for the Sedov blast example..."
python makeIC.py
fi
../swift -s -t 16 sedov.yml
# Run SWIFT
../swift -s -t 4 sedov.yml
# Plot the solution
python plotSolution.py 5
examples/SedovBlast_3D/sedov.py deleted 100755 → 0
View file @ d5f913b8
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
#
# 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 numpy as np
import scipy.integrate as integrate
import scipy.optimize as optimize
import os
# Calculate the analytical solution of the Sedov-Taylor shock wave for a given
# number of dimensions, gamma and time. We assume dimensionless units and a
# setup with constant density 1, pressure and velocity 0. An energy 1 is
# inserted in the center at t 0.
#
# The solution is a self-similar shock wave, which was described in detail by
# Sedov (1959). We follow his notations and solution method.
#
# The position of the shock at time t is given by
# r2 = (E/rho1)^(1/(2+nu)) * t^(2/(2+nu))
# the energy E is related to the inserted energy E0 by E0 = alpha*E, with alpha
# a constant which has to be calculated by imposing energy conservation.
#
# The density for a given radius at a certain time is determined by introducing
# the dimensionless position coordinate lambda = r/r2. The density profile as
# a function of lambda is constant in time and given by
# rho = rho1 * R(V, gamma, nu)
# and
# V = V(lambda, gamma, nu)
#
# The function V(lambda, gamma, nu) is found by solving a differential equation
# described in detail in Sedov (1959) chapter 4, section 5. Alpha is calculated
# from the integrals in section 11 of the same chapter.
#
# Numerically, the complete solution requires the use of 3 quadratures and 1
# root solver, which are implemented using the GNU Scientific Library (GSL).
# Since some quadratures call functions that themselves contain quadratures,
# the problem is highly non-linear and complex and takes a while to converge.
# Therefore, we tabulate the alpha values and profile values the first time
# a given set of gamma and nu values is requested and reuse these tabulated
# values.
#
# Reference:
# Sedov (1959): Sedov, L., Similitude et dimensions en mecanique (7th ed.;
# Moscou: Editions Mir) - french translation of the original
# book from 1959.
# dimensionless variable z = gamma*P/R as a function of V (and gamma and nu)
# R is a dimensionless density, while P is a dimensionless pressure
# z is hence a sort of dimensionless sound speed
# The expression below corresponds to eq. 11.9 in Sedov (1959), chapter 4
def _z(V, gamma, nu):
if V == 2./(nu+2.)/gamma:
return 0.
else:
return (gamma-1.)*V*V*(V-2./(2.+nu))/2./(2./(2.+nu)/gamma-V)
# differential equation that needs to be solved to obtain lambda for a given V
# corresponds to eq. 5.11 in Sedov (1959), chapter 4 (omega = 0)
def _dlnlambda_dV(V, gamma, nu):
nom = _z(V, gamma, nu) - (V-2./(nu+2.))*(V-2./(nu+2.))
denom = V*(V-1.)*(V-2./(nu+2.))+nu*(2./(nu+2.)/gamma-V)*_z(V, gamma, nu)
return nom/denom
# dimensionless variable lambda = r/r2 as a function of V (and gamma and nu)
# found by solving differential equation 5.11 in Sedov (1959), chapter 4
# (omega = 0)
def _lambda(V, gamma, nu):
if V == 2./(nu+2.)/gamma:
return 0.
else:
V0 = 4./(nu+2.)/(gamma+1.)
integral, err = integrate.quad(_dlnlambda_dV, V, V0, (gamma, nu),
limit = 8000)
return np.exp(-integral)
# dimensionless variable R = rho/rho1 as a function of V (and gamma and nu)
# found by inverting eq. 5.12 in Sedov (1959), chapter 4 (omega = 0)
# the integration constant C is found by inserting the R, V and z values
# at the shock wave, where lambda is 1. These correspond to eq. 11.8 in Sedov
# (1959), chapter 4.
def _R(V, gamma, nu):
if V == 2./(nu+2.)/gamma:
return 0.
else:
C = 8.*gamma*(gamma-1.)/(nu+2.)/(nu+2.)/(gamma+1.)/(gamma+1.) \
*((gamma-1.)/(gamma+1.))**(gamma-2.) \
*(4./(nu+2.)/(gamma+1.)-2./(nu+2.))
lambda1 = _lambda(V, gamma, nu)
lambda5 = lambda1**(nu+2)
return (_z(V, gamma, nu)*(V-2./(nu+2.))*lambda5/C)**(1./(gamma-2.))
# function of which we need to find the zero point to invert lambda(V)
def _lambda_min_lambda(V, lambdax, gamma, nu):
return _lambda(V, gamma, nu) - lambdax
# dimensionless variable V = v*t/r as a function of lambda (and gamma and nu)
# found by inverting the function lamdba(V) which is found by solving
# differential equation 5.11 in Sedov (1959), chapter 4 (omega = 0)
# the invertion is done by searching the zero point of the function
# lambda_min_lambda defined above
def _V_inv(lambdax, gamma, nu):
if lambdax == 0.:
return 2./(2.+nu)/gamma;
else:
return optimize.brentq(_lambda_min_lambda, 2./(nu+2.)/gamma,
4./(nu+2.)/(gamma+1.), (lambdax, gamma, nu))
# integrand of the first integral in eq. 11.24 in Sedov (1959), chapter 4
def _integrandum1(lambdax, gamma, nu):
V = _V_inv(lambdax, gamma, nu)
if nu == 2:
return _R(V, gamma, nu)*V**2*lambdax**3
else:
return _R(V, gamma, nu)*V**2*lambdax**4
# integrand of the second integral in eq. 11.24 in Sedov (1959), chapter 4
def _integrandum2(lambdax, gamma, nu):
V = _V_inv(lambdax, gamma, nu)
if V == 2./(nu+2.)/gamma:
P = 0.
else:
P = _z(V, gamma, nu)*_R(V, gamma, nu)/gamma
if nu == 2:
return P*lambdax**3
else:
return P*lambdax**4
# calculate alpha = E0/E
# this corresponds to eq. 11.24 in Sedov (1959), chapter 4
def get_alpha(gamma, nu):
integral1, err1 = integrate.quad(_integrandum1, 0., 1., (gamma, nu))
integral2, err2 = integrate.quad(_integrandum2, 0., 1., (gamma, nu))