diff --git a/examples/SedovBlast_2D/makeIC.py b/examples/SedovBlast_2D/makeIC.py
index a30e727b4dd42c1f058827959cf12e3b4f152181..05233576f63f90b8d448aaa75fa6bfe7fce1f0e8 100644
--- a/examples/SedovBlast_2D/makeIC.py
+++ b/examples/SedovBlast_2D/makeIC.py
@@ -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
diff --git a/examples/SedovBlast_3D/makeIC.py b/examples/SedovBlast_3D/makeIC.py
index e64942e8e92ee6fe67142f841f566019b1a668be..3c1e36a74b53ece8b886e2fcbe5d9178d9deefbc 100644
--- a/examples/SedovBlast_3D/makeIC.py
+++ b/examples/SedovBlast_3D/makeIC.py
@@ -1,7 +1,6 @@
###############################################################################
# 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')
diff --git a/examples/SedovBlast_3D/makeIC_fcc.py b/examples/SedovBlast_3D/makeIC_fcc.py
deleted file mode 100644
index 0d3a017a9b7f3b30b61e723e3d1646d7797b40a4..0000000000000000000000000000000000000000
--- a/examples/SedovBlast_3D/makeIC_fcc.py
+++ /dev/null
@@ -1,122 +0,0 @@
-###############################################################################
- # 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()
diff --git a/examples/SedovBlast_3D/profile.py b/examples/SedovBlast_3D/profile.py
deleted file mode 100644
index 1b1dfb389b2c9f6d0e5af02a17a07068f03bfd92..0000000000000000000000000000000000000000
--- a/examples/SedovBlast_3D/profile.py
+++ /dev/null
@@ -1,46 +0,0 @@
-###############################################################################
- # 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()
diff --git a/examples/SedovBlast_3D/rdf.py b/examples/SedovBlast_3D/rdf.py
deleted file mode 100644
index 7f932cc814dc36e14e1bef52e33cf5ed1f527dfd..0000000000000000000000000000000000000000
--- a/examples/SedovBlast_3D/rdf.py
+++ /dev/null
@@ -1,121 +0,0 @@
-###############################################################################
- # 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] )
-
-
diff --git a/examples/SedovBlast_3D/run.sh b/examples/SedovBlast_3D/run.sh
index f71830eb6a66a9ea84e93fd1bed1261b1cb42b7b..afc9ff257eff334d700dbd139c7ffbf0225fbf7b 100755
--- a/examples/SedovBlast_3D/run.sh
+++ b/examples/SedovBlast_3D/run.sh
@@ -1,10 +1,19 @@
#!/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
diff --git a/examples/SedovBlast_3D/sedov.py b/examples/SedovBlast_3D/sedov.py
deleted file mode 100755
index ca908fe65cd15ddd2b0e4fb41a402492432a0690..0000000000000000000000000000000000000000
--- a/examples/SedovBlast_3D/sedov.py
+++ /dev/null
@@ -1,212 +0,0 @@
-###############################################################################
- # 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))
-
- if nu == 2:
- return np.pi*integral1+2.*np.pi/(gamma-1.)*integral2
- else:
- return 2.*np.pi*integral1+4.*np.pi/(gamma-1.)*integral2
-
-# get the analytical solution for the Sedov-Taylor blastwave given an input
-# energy E, adiabatic index gamma, and number of dimensions nu, at time t and
-# with a maximal outer region radius maxr
-def get_analytical_solution(E, gamma, nu, t, maxr = 1.):
- # we check for the existance of a datafile with precomputed alpha and
- # profile values
- # if it does not exist, we calculate it here and write it out
- # calculation of alpha and the profile takes a while...
- lvec = np.zeros(1000)
- Rvec = np.zeros(1000)
- fname = "sedov_profile_gamma_{gamma}_nu_{nu}.dat".format(gamma = gamma,
- nu = nu)
- if os.path.exists(fname):
- file = open(fname, "r")
- lines = file.readlines()
- alpha = float(lines[0])
- for i in range(len(lines)-1):
- data = lines[i+1].split()
- lvec[i] = float(data[0])
- Rvec[i] = float(data[1])
- else:
- alpha = get_alpha(gamma, nu)
- for i in range(1000):
- lvec[i] = (i+1)*0.001
- V = _V_inv(lvec[i], gamma, nu)
- Rvec[i] = _R(V, gamma, nu)
- file = open(fname, "w")
- file.write("{alpha}\n".format(alpha = alpha))
- for i in range(1000):
- file.write("{l}\t{R}\n".format(l = lvec[i], R = Rvec[i]))
-
- xvec = np.zeros(1002)
- rhovec = np.zeros(1002)
- if nu == 2:
- r2 = (E/alpha)**0.25*np.sqrt(t)
- else:
- r2 = (E/alpha)**0.2*t**0.4
-
- for i in range(1000):
- xvec[i] = lvec[i]*r2
- rhovec[i] = Rvec[i]
- xvec[1000] = 1.001*r2
- xvec[1001] = maxr
- rhovec[1000] = 1.
- rhovec[1001] = 1.
-
- return xvec, rhovec
-
-def main():
- E = 1.
- gamma = 1.66667
- nu = 3
- t = 0.001
- x, rho = get_analytical_solution(E, gamma, nu, t)
- for i in range(len(x)):
- print x[i], rho[i]
-
-if __name__ == "__main__":
- main()
diff --git a/examples/SedovBlast_3D/sedov.yml b/examples/SedovBlast_3D/sedov.yml
index eb95fd85d5599145b2ff037256dcde72fc306209..6f519835d26ff5aa851ffb8999e650815c522cd3 100644
--- a/examples/SedovBlast_3D/sedov.yml
+++ b/examples/SedovBlast_3D/sedov.yml
@@ -9,15 +9,15 @@ 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: 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).
- dt_max: 1e-2 # The maximal 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: sedov # 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)
+ delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
@@ -27,7 +27,7 @@ Statistics:
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: 1. # Maximal smoothing length allowed (in internal units).
+ 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
diff --git a/examples/SedovBlast_3D/solution.py b/examples/SedovBlast_3D/solution.py
deleted file mode 100644
index 9335e22bdd76585e8d53d3dba9f9a435197f92fc..0000000000000000000000000000000000000000
--- a/examples/SedovBlast_3D/solution.py
+++ /dev/null
@@ -1,114 +0,0 @@
-###############################################################################
- # 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 3D Sedov blast wave.
-# The script works for a given initial box and dumped energy 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 radial points between r=0 and r=R_max.
-# Follows the solution in Landau & Lifschitz
-
-# Parameters
-rho_0 = 1. # Background Density
-P_0 = 1.e-5 # Background Pressure
-E_0 = 1.e2 # Energy of the explosion
-gamma = 5./3. # Gas polytropic index
-
-t = 0.275 # Time of the solution
-
-N = 1000 # Number of radial points
-R_max = 3. # Maximal radius
-
-# ---------------------------------------------------------------
-# Don't touch anything after this.
-# ---------------------------------------------------------------
-
-if gamma == 5./3.:
- alpha = 0.49
-else:
- print "Unknown value for alpha"
- exit
-
-# Position and velocity of the shock
-r_shock = (E_0 / (alpha * rho_0))**(1./5.) * t**(2./5.)
-v_shock = (1./5.) * (1./alpha)**(1./5.) * ((E_0 * t**2. / rho_0)**(-4./5.)) * E_0 * (t / rho_0)
-
-# Prepare arrays
-delta_r = R_max / N
-r_s = arange(0, R_max, delta_r)
-rho_s = ones(N) * rho_0
-P_s = ones(N) * P_0
-u_s = ones(N)
-v_s = zeros(N)
-
-# State on the shock
-rho_shock = rho_0 * (gamma+1.)/(gamma-1.)
-P_shock = 2./(gamma+1.) * rho_shock * v_shock**2
-
-# Integer position of the shock
-i_shock = min(floor(r_shock /delta_r), N)
-
-# Dimensionless velocity and its spatial derivative
-v_bar0 = (gamma+1.)/(2.*gamma)
-deltaV_bar = (1.0 - v_bar0) / (i_shock - 1)
-
-def rho_dimensionless(v_bar):
- power1 = (2./(gamma-2.))
- power2 = -(12.-7.*gamma+13.*gamma**2)/(2.-3.*gamma-11.*gamma**2+6.*gamma**3)
- power3 = 3./(1.+2.*gamma)
- term1 = ((1.+ gamma - 2.*v_bar)/(gamma-1.))**power1
- term2 = ((5.+5.*gamma+2.*v_bar-6.*gamma*v_bar)/(7.-gamma))**power2
- term3 = ((2.*gamma*v_bar - gamma -1.)/(gamma-1.))**power3
- return term1 * term2 * term3
-
-def P_dimensionless(v_bar):
- return (gamma+1. - 2.*v_bar)/(2.*gamma*v_bar - gamma - 1.)*v_bar**2 * rho_dimensionless(v_bar)
-
-def r_dimensionless(v_bar):
- power1 = (-12.+7.*gamma-13.*gamma**2)/(5.*(-1.+gamma+6.*gamma**2))
- power2 = (gamma - 1.)/(1.+2.*gamma)
- term1 = ((5.+5.*gamma+2.*v_bar-6.*gamma*v_bar)/(7.-gamma))**power1
- term2 = ((2.*gamma*v_bar-gamma-1.)/(gamma-1.))**power2
- return v_bar**(-2./5.)*term1*term2
-
-# Generate solution
-for i in range(1,int(i_shock)):
- v_bar = v_bar0 + (i-1)*deltaV_bar
- r_s[i] = r_dimensionless(v_bar) * r_shock
- rho_s[i] = rho_dimensionless(v_bar) * rho_shock
- P_s[i] = P_shock * (r_s[i]/r_shock)**2 * P_dimensionless(v_bar)
- v_s[i] = (4./(5.*(gamma+1.)) * r_s[i] / t * v_bar)
-
-u_s = P_s / ((gamma - 1.)*rho_s) # internal energy
-s_s = P_s / rho_s**gamma # entropic function
-rho_s[0] = 0.
-P_s[0] = P_s[1] # dirty...
-
-
-savetxt("rho.dat", column_stack((r_s, rho_s)))
-savetxt("P.dat", column_stack((r_s, P_s)))
-savetxt("v.dat", column_stack((r_s, v_s)))
-savetxt("u.dat", column_stack((r_s, u_s)))
-savetxt("s.dat", column_stack((r_s, s_s)))
-
-
-
diff --git a/examples/SodShock_3D/makeIC.py b/examples/SodShock_3D/makeIC.py
index 8ae19050c11c0712579b44646c8870d7574d113b..28c2cfd82640c9d3a040d8d58ba31154c0719075 100644
--- a/examples/SodShock_3D/makeIC.py
+++ b/examples/SodShock_3D/makeIC.py
@@ -19,7 +19,6 @@
##############################################################################
import h5py
-import random
from numpy import *
# Generates a swift IC file for the Sod Shock in a periodic box
diff --git a/examples/SodShock_3D/sodShock.yml b/examples/SodShock_3D/sodShock.yml
index 19bacb57f3bfb173063dbac5d752929763dede29..a46d521511e9885ba2e5425ce1aa730403be533a 100644
--- a/examples/SodShock_3D/sodShock.yml
+++ b/examples/SodShock_3D/sodShock.yml
@@ -17,7 +17,7 @@ TimeIntegration:
Snapshots:
basename: sod # 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)
+ delta_time: 0.05 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics: