Skip to content
Snippets Groups Projects
Commit b2345b3c authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' into disk-patch

parents 911be52a 954b572d
No related branches found
No related tags found
1 merge request!217Disk patch
Showing
with 808 additions and 504 deletions
......@@ -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
......
......@@ -469,7 +469,7 @@ if test "$enable_warn" != "no"; then
CFLAGS="$CFLAGS -Wall"
;;
intel)
CFLAGS="$CFLAGS -w2"
CFLAGS="$CFLAGS -w2 -Wunused-variable"
;;
*)
AX_CFLAGS_WARN_ALL
......
###############################################################################
# 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)))
#!/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.02 # 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
###############################################################################
# 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
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 = "glassPlane_128.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**2 / 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, 0.2]
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]
#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.
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 = 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.reshape((numPart,1))
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h.reshape((numPart,1))
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u.reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
ds[()] = ids.reshape((numPart,1))
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 glassPlane_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
# 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.5 # 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: kelvinHelmholtz # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.25 # 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.01 # 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: ./kelvinHelmholtz.hdf5 # The file to read
###############################################################################
# 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
from numpy import *
import sys
# Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
# Parameters
L2 = 128 # Particles along one edge in the low-density region
gamma = 5./3. # Gas adiabatic index
P1 = 2.5 # Central region pressure
P2 = 2.5 # Outskirts pressure
v1 = 0.5 # Central region velocity
v2 = -0.5 # Outskirts vlocity
rho1 = 2 # Central density
rho2 = 1 # Outskirts density
omega0 = 0.1
sigma = 0.05 / sqrt(2)
fileOutputName = "kelvinHelmholtz.hdf5"
#---------------------------------------------------
# Start by generating grids of particles at the two densities
numPart2 = L2 * L2
L1 = int(sqrt(numPart2 / rho2 * rho1))
numPart1 = L1 * L1
#print "N2 =", numPart2, "N1 =", numPart1
#print "L2 =", L2, "L1 = ", L1
#print "rho2 =", rho2, "rho1 =", (float(L1*L1)) / (float(L2*L2))
coords1 = zeros((numPart1, 3))
coords2 = zeros((numPart2, 3))
h1 = ones(numPart1) * 1.2348 / L1
h2 = ones(numPart2) * 1.2348 / L2
m1 = zeros(numPart1)
m2 = zeros(numPart2)
u1 = zeros(numPart1)
u2 = zeros(numPart2)
vel1 = zeros((numPart1, 3))
vel2 = zeros((numPart2, 3))
# Particles in the central region
for i in range(L1):
for j in range(L1):
index = i * L1 + j
x = i / float(L1) + 1. / (2. * L1)
y = j / float(L1) + 1. / (2. * L1)
coords1[index, 0] = x
coords1[index, 1] = y
u1[index] = P1 / (rho1 * (gamma-1.))
vel1[index, 0] = v1
# Particles in the outskirts
for i in range(L2):
for j in range(L2):
index = i * L2 + j
x = i / float(L2) + 1. / (2. * L2)
y = j / float(L2) + 1. / (2. * L2)
coords2[index, 0] = x
coords2[index, 1] = y
u2[index] = P2 / (rho2 * (gamma-1.))
vel2[index, 0] = v2
# Now concatenate arrays
where1 = abs(coords1[:,1]-0.5) < 0.25
where2 = abs(coords2[:,1]-0.5) > 0.25
coords = append(coords1[where1, :], coords2[where2, :], axis=0)
#print L2*(L2/2), L1*(L1/2)
#print shape(coords), shape(coords1[where1,:]), shape(coords2[where2,:])
#print shape(coords), shape(logical_not(coords1[where1,:])), shape(logical_not(coords2[where2,:]))
vel = append(vel1[where1, :], vel2[where2, :], axis=0)
h = append(h1[where1], h2[where2], axis=0)
m = append(m1[where1], m2[where2], axis=0)
u = append(u1[where1], u2[where2], axis=0)
numPart = size(h)
ids = linspace(1, numPart, numPart)
m[:] = (0.5 * rho1 + 0.5 * rho2) / float(numPart)
# Velocity perturbation
vel[:,1] = omega0 * sin(4*pi*coords[:,0]) * (exp(-(coords[:,1]-0.25)**2 / (2 * sigma**2)) + exp(-(coords[:,1]-0.75)**2 / (2 * sigma**2)))
#File
fileOutput = h5py.File(fileOutputName, 'w')
# Header
grp = fileOutput.create_group("/Header")
grp.attrs["BoxSize"] = [1., 1., 0.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]
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]
#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.
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 = fileOutput.create_group("/PartType0")
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
ds[()] = coords
ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = vel
ds = grp.create_dataset('Masses', (numPart, 1), 'f')
ds[()] = m.reshape((numPart,1))
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h.reshape((numPart,1))
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u.reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
ds[()] = ids.reshape((numPart,1))
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
P1 = 2.5 # Central region pressure
P2 = 2.5 # Outskirts pressure
v1 = 0.5 # Central region velocity
v2 = -0.5 # Outskirts vlocity
rho1 = 2 # Central density
rho2 = 1 # Outskirts density
# ---------------------------------------------------------------
# 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("kelvinHelmholtz_%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"][:,:]
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)
scatter(pos[:,0], pos[:,1], c=vel[:,0], cmap="PuBu", edgecolors='face', s=4, vmin=-1., vmax=1.)
text(0.97, 0.97, "${\\rm{Velocity~along}}~x$", ha="right", va="top", backgroundcolor="w")
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Position}}~y$", labelpad=0)
xlim(0, 1)
ylim(0, 1)
# Radial density profile --------------------------------
subplot(232)
scatter(pos[:,0], pos[:,1], c=rho, cmap="PuBu", edgecolors='face', s=4, vmin=0.8, vmax=2.2)
text(0.97, 0.97, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Position}}~y$", labelpad=0)
xlim(0, 1)
ylim(0, 1)
# Radial pressure profile --------------------------------
subplot(233)
scatter(pos[:,0], pos[:,1], c=P, cmap="PuBu", edgecolors='face', s=4, vmin=1, vmax=4)
text(0.97, 0.97, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Position}}~y$", labelpad=0)
xlim(0, 1)
ylim(0, 1)
# Internal energy profile --------------------------------
subplot(234)
scatter(pos[:,0], pos[:,1], c=u, cmap="PuBu", edgecolors='face', s=4, vmin=1.5, vmax=5.)
text(0.97, 0.97, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w")
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Position}}~y$", labelpad=0)
xlim(0, 1)
ylim(0, 1)
# Radial entropy profile --------------------------------
subplot(235)
scatter(pos[:,0], pos[:,1], c=S, cmap="PuBu", edgecolors='face', s=4, vmin=0.5, vmax=3.)
text(0.97, 0.97, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Position}}~y$", labelpad=0)
xlim(0, 1)
ylim(0, 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, "Kelvin-Helmholtz instability at $t=%.2f$"%(time), fontsize=10)
text(-0.49, 0.8, "Centre:~~~ $(P, \\rho, v) = (%.3f, %.3f, %.3f)$"%(P1, rho1, v1), fontsize=10)
text(-0.49, 0.7, "Outskirts: $(P, \\rho, v) = (%.3f, %.3f, %.3f)$"%(P2, rho2, v2), 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("KelvinHelmholtz.png", dpi=200)
#!/bin/bash
# Generate the initial conditions if they are not present.
if [ ! -e kelvinHelmholtz.hdf5 ]
then
echo "Generating initial conditions for the Kelvin-Helmholtz example..."
python makeIC.py
fi
# Run SWIFT
../swift -s -t 1 kelvinHelmholtz.yml
# Plot the solution
python plotSolution.py 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,28 +18,28 @@
##############################################################################
import h5py
import sys
import random
from numpy import *
# Generates a swift IC file for the Gresho Vortex in a periodic box
# 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
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
P0 = 0. # Constant additional pressure (should have no impact on the dynamics)
fileName = "greshoVortex.hdf5"
vol = boxSize[0] * boxSize[1] * boxSize[2]
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**3 / factor
mass = boxSize[0]*boxSize[1]*boxSize[2] * rho / numPart
numPart = L**2
mass = boxSize**2 * rho / numPart
internalEnergy = P / ((gamma - 1.)*rho)
#Generate particles
coords = zeros((numPart, 3))
......@@ -50,43 +49,26 @@ 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
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')
......@@ -101,6 +83,7 @@ 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, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total"] = numPart
#Runtime parameters
grp = file.create_group("/RuntimePars")
......@@ -126,9 +109,7 @@ ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
ds[()] = ids[:]
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
ds[()] = ids + 1
file.close()
# 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: 10. # 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: perturbedPlane # 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-3 # 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: ./perturbedPlane.hdf5 # The file to read
File moved
###############################################################################
# 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()
###############################################################################
# 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] )
###############################################################################
# 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()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment