diff --git a/examples/SedovBlast_2D/makeIC.py b/examples/SedovBlast_2D/makeIC.py index dc2c574bb76387ef3bebd714423e4612e5838e0d..a30e727b4dd42c1f058827959cf12e3b4f152181 100644 --- a/examples/SedovBlast_2D/makeIC.py +++ b/examples/SedovBlast_2D/makeIC.py @@ -31,8 +31,6 @@ E0= 1. # Energy of the explosion N_inject = 21 # Number of particles in which to inject energy fileName = "sedov.hdf5" -#L = 101 - #--------------------------------------------------- glass = h5py.File("glassPlane_128.hdf5", "r") @@ -50,10 +48,9 @@ m = zeros(numPart) u = zeros(numPart) r = zeros(numPart) -for i in range(numPart): - r[i] = sqrt((pos[i,0] - 0.5)**2 + (pos[i,1] - 0.5)**2) - m[i] = rho0 * vol / numPart - u[i] = P0 / (rho0 * (gamma - 1)) +r = sqrt((pos[:,0] - 0.5)**2 + (pos[:,1] - 0.5)**2) +m[:] = rho0 * vol / numPart +u[:] = P0 / (rho0 * (gamma - 1)) # Make the central particle detonate index = argsort(r) diff --git a/examples/SedovBlast_2D/plotSolution.py b/examples/SedovBlast_2D/plotSolution.py index af7c3b525dc25b2b10399f9cf27553999b28ab92..69e4e1232dd5c14b06e8a705f4add391f1f597f0 100644 --- a/examples/SedovBlast_2D/plotSolution.py +++ b/examples/SedovBlast_2D/plotSolution.py @@ -18,7 +18,7 @@ # ############################################################################## -# Computes the analytical solution of the 3D Sedov blast wave. +# Computes the analytical solution of the 2D Sedov blast wave. # The script works for a given initial box and dumped energy and computes the solution at a later time t. # Parameters @@ -85,6 +85,8 @@ P = sim["/PartType0/Pressure"][:] rho = sim["/PartType0/Density"][:] +# Now, work our the solution.... + from scipy.special import gamma as Gamma from numpy import * @@ -190,6 +192,8 @@ def sedov(t, E0, rho0, g, n=1000, nu=3): rho *= rho0 return r, p, rho, u, r_s, p_s, rho_s, u_s, shock_speed + +# The main properties of the solution r_s, P_s, rho_s, v_s, r_shock, _, _, _, _ = sedov(time, E_0, rho_0, gas_gamma, 1000, 2) # Append points for after the shock @@ -202,6 +206,8 @@ v_s = np.insert(v_s, np.size(v_s), [0, 0]) u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy s_s = P_s / rho_s**gas_gamma # entropic function + + # Plot the interesting quantities figure() diff --git a/examples/SedovBlast_2D/sedov.py b/examples/SedovBlast_2D/sedov.py deleted file mode 100755 index 2439a7fda91831a7c54350a597db5f026e36fee9..0000000000000000000000000000000000000000 --- a/examples/SedovBlast_2D/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 = 2 - 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()