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

Updated the script generating the solution of the Sedov blast. It now outputs... 

Updated the script generating the solution of the Sedov blast. It now outputs all the relevant arrays. The script as is corresponds to the ICs in makeIC.py.


Former-commit-id: bc3e5383ca0c4ffdbdc7368ca153b96cd8be155d
parent c773c6cc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
examples/SedovBlast/solution.py
+ 114
− 214
View file @ 16cbcb7b
"""
Peter Creasey
p.e.creasey.00@googlemail.com
solution to the Sedov problem
based on the C code by Aamer Haque
"""
from scipy.special import gamma as Gamma
from numpy import power, arange, empty, float64, log, exp, pi, diff, inner, outer, array
def calc_a(g,nu=3):
"""
exponents of the polynomials of the sedov solution
g - the polytropic gamma
nu - the dimension
"""
a = [0]*8
a[0] = 2.0 / (nu + 2)
a[2] = (1-g) / (2*(g-1) + nu)
a[3] = nu / (2*(g-1) + nu)
a[5] = 2 / (g-2)
a[6] = g / (2*(g-1) + nu)
a[1] = (((nu+2)*g)/(2.0+nu*(g-1.0)) ) * ( (2.0*nu*(2.0-g))/(g*(nu+2.0)**2) - a[2])
a[4] = a[1]*(nu+2) / (2-g)
a[7] = (2 + nu*(g-1))*a[1]/(nu*(2-g))
return a
def calc_beta(v, g, nu=3):
"""
beta values for the sedov solution (coefficients of the polynomials of the similarity variables)
v - the similarity variable
g - the polytropic gamma
nu- the dimension
"""
beta = (nu+2) * (g+1) * array((0.25, (g/(g-1))*0.5,
-(2 + nu*(g-1))/2.0 / ((nu+2)*(g+1) -2*(2 + nu*(g-1))),
-0.5/(g-1)), dtype=float64)
beta = outer(beta, v)
beta += (g+1) * array((0.0, -1.0/(g-1),
(nu+2) / ((nu+2)*(g+1) -2.0*(2 + nu*(g-1))),
1.0/(g-1)), dtype=float64).reshape((4,1))
return beta
def sedov(t, E0, rho0, g, n=1000, nu=3):
"""
solve the sedov problem
t - the time
E0 - the initial energy
rho0 - the initial density
n - number of points (10000)
nu - the dimension
g - the polytropic gas gamma
"""
# the similarity variable
v_min = 2.0 / ((nu + 2) * g)
v_max = 4.0 / ((nu + 2) * (g + 1))
v = v_min + arange(n) * (v_max - v_min) / (n - 1.0)
a = calc_a(g, nu)
beta = calc_beta(v, g=g, nu=nu)
lbeta = log(beta)
r = exp(-a[0] * lbeta[0] - a[2] * lbeta[1] - a[1] * lbeta[2])
rho = ((g + 1.0) / (g - 1.0)) * exp(a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2])
p = exp(nu * a[0] * lbeta[0] + (a[5] + 1) * lbeta[3] + (a[4] - 2 * a[1]) * lbeta[2])
u = beta[0] * r * 4.0 / ((g + 1) * (nu + 2))
p *= 8.0 / ((g + 1) * (nu + 2) * (nu + 2))
# we have to take extra care at v=v_min, since this can be a special point.
# It is not a singularity, however, the gradients of our variables (wrt v) are.
# r -> 0, u -> 0, rho -> 0, p-> constant
u[0] = 0.0; rho[0] = 0.0; r[0] = 0.0; p[0] = p[1]
# volume of an n-sphere
vol = (pi ** (nu / 2.0) / Gamma(nu / 2.0 + 1)) * power(r, nu)
# note we choose to evaluate the integral in this way because the
# volumes of the first few elements (i.e near v=vmin) are shrinking
# very slowly, so we dramatically improve the error convergence by
# finding the volumes exactly. This is most important for the
# pressure integral, as this is on the order of the volume.
# (dimensionless) energy of the model solution
de = rho * u * u * 0.5 + p / (g - 1)
# integrate (trapezium rule)
q = inner(de[1:] + de[:-1], diff(vol)) * 0.5
# the factor to convert to this particular problem
fac = (q * (t ** nu) * rho0 / E0) ** (-1.0 / (nu + 2))
# shock speed
shock_speed = fac * (2.0 / (nu + 2))
rho_s = ((g + 1) / (g - 1)) * rho0
r_s = shock_speed * t * (nu + 2) / 2.0
p_s = (2.0 * rho0 * shock_speed * shock_speed) / (g + 1)
u_s = (2.0 * shock_speed) / (g + 1)
r *= fac * t
u *= fac
p *= fac * fac * rho0
rho *= rho0
return r, p, rho, u, r_s, p_s, rho_s, u_s, shock_speed
def test():
""" draw a 3d sedov solution """
import pylab as pl
gamma = 5.0/3.0
r,p,rho,u,r_s,p_s,rho_s,u_s,shock_speed = \
sedov(t=0.05, E0=5.0, rho0=5.0, g=gamma)
print 'rho shock', rho_s
print 'p shock', p_s
print 'u shock', u_s
print 'r shock', r_s
print 'Dimensionless var (E/rho) t^2 r^-5', (5.0 /5.0)* 0.05**0.4 * r[-1]**-1.0
vols = (4/3.0)*pi*r*r*r
dv = vols.copy()
dv[1:] = diff(dv)
# thermal and kinetic energy
te = (p*dv/(gamma-1))
ke = (rho*u*u*0.5*dv)
energy = te.sum() + ke.sum()
mass = 0.5*inner(rho[1:]+rho[:-1],dv[1:])
print 'density', mass / (4/3.0 * pi * r_s**3)
print 'energy', energy
print 'shock speed', shock_speed
pl.plot(r/r_s,rho/rho_s, label=r'$\rho/\rho_s$')
pl.plot(r/r_s,p/p_s, label=r'$p/p_s$')
pl.plot(r/r_s,u/u_s, label=r'$u/u_s$')
pl.legend(loc='upper left')
pl.show()
def test2():
""" test momentum and mass conservation in 3d """
import pylab as pl
r,p,rho,u,r_s,p_s,rho_s,u_s,shock_speed = \
sedov(t=0.05, E0=5.0, rho0=5.0, g=5.0/3.0,n=10000)
dt = 1e-5
r2,p2,rho2,u2 = sedov(t=0.05+dt, E0=5.0, rho0=5.0, g=5.0/3.0, n=9000)[:4]
# align the results
from numpy import interp, gradient
p2 = interp(r,r2,p2)
rho2 = interp(r,r2,rho2)
u2 = interp(r,r2,u2)
# mass conservation
pl.plot(r, -gradient(rho*u*r*r)/(r*r*gradient(r)), 'b', label=r'$\frac{1}{r^2}\frac{\partial}{\partial r} \rho u r^2$')
pl.plot(r, (rho2-rho)/dt, 'k', label=r'$\frac{\partial \rho}{\partial t}$')
# momentum conservation
pl.plot(r, -gradient(p)/gradient(r), 'g',label=r'$-\frac{\partial p}{\partial r}$')
pl.plot(r, rho*((u2-u)/dt+u*gradient(u)/gradient(r)), 'r',label=r'$\rho \left( \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial r} \right)$')
pl.legend(loc='lower left')
pl.show()
def test3():
""" draw a 2d sedov solution """
import pylab as pl
r,p,rho,u,r_s,p_s,rho_s,u_s,shock_speed = sedov(t=1.2, E0=1, rho0=1, g=5.0/3.0, nu=2)
print 'rho shock', rho_s
print 'p shock', p_s
print 'u shock', u_s
print 'r shock', r_s
area = pi*r*r
dv = area.copy()
dv[1:] = diff(dv)
# thermal and kinetic energy
te = (p*dv/(5.0/3.0-1))
ke = (rho*u*u*0.5*dv)
#pl.plot(arange(te.size), ke, 'x')
#pl.show()
print 'r0', r[:2]
energy = te.sum() + ke.sum()
mass = 0.5*inner(rho[1:]+rho[:-1],dv[1:])
print 'density', mass / (pi * r_s**2)
print 'energy', energy
print 'shock speed', shock_speed
#pl.plot(r/r_s,rho/rho_s, 'b,',label=r'$\rho/\rho_s$')
#pl.plot(r/r_s,p/p_s,'r',label=r'$p/p_s$')
#pl.plot(r/r_s,u/u_s, 'g,',label=r'$u/u_s$')
# pl.plot(r,rho,'b',label='$\\rho$')
pl.plot(r,u,'b',label='$u$')
# pl.plot(r,p,'b',label='$P$')
pl.legend(loc='upper left')
# pl.xlim(0,2)
# pl.ylim(0,5)
pl.show()
if __name__=='__main__':
test()
#test2()
# test3()
###############################################################################
# This file is part of SWIFT.
# Coypright (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.e5 # Energy of the explosion
gamma = 5./3. # Gas polytropic index
t = 0.1 # Time of the solution
N = 1000 # Number of radial points
R_max = 5. # 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)))
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment