Self-similar Sedov-Taylor blast wave solution

sedov-taylor.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Feb  5 13:54:03 2025

@author: alankar
"""

import numpy as np
from scipy.integrate import solve_ivp, simpson
import matplotlib.pyplot as plt
import matplotlib

dark = True

## Plot Styling
matplotlib.rcParams["xtick.direction"] = "in"
matplotlib.rcParams["ytick.direction"] = "in"
matplotlib.rcParams["xtick.top"] = False
matplotlib.rcParams["ytick.right"] = True
matplotlib.rcParams["xtick.minor.visible"] = True
matplotlib.rcParams["ytick.minor.visible"] = True
matplotlib.rcParams["axes.grid"] = True
matplotlib.rcParams["grid.linestyle"] = ":"
matplotlib.rcParams["grid.linewidth"] = 1.0
matplotlib.rcParams["grid.color"] = "gray" if not(dark) else "white"
matplotlib.rcParams["grid.alpha"] = 0.3 if not(dark) else 0.8
matplotlib.rcParams["lines.dash_capstyle"] = "round"
matplotlib.rcParams["lines.solid_capstyle"] = "round"
matplotlib.rcParams["legend.handletextpad"] = 0.4
matplotlib.rcParams["axes.linewidth"] = 1.0
matplotlib.rcParams["lines.linewidth"] = 3.5
matplotlib.rcParams["ytick.major.width"] = 1.2
matplotlib.rcParams["xtick.major.width"] = 1.2
matplotlib.rcParams["ytick.minor.width"] = 1.0
matplotlib.rcParams["xtick.minor.width"] = 1.0
matplotlib.rcParams["ytick.major.size"] = 11.0
matplotlib.rcParams["xtick.major.size"] = 11.0
matplotlib.rcParams["ytick.minor.size"] = 5.0
matplotlib.rcParams["xtick.minor.size"] = 5.0
matplotlib.rcParams["xtick.major.pad"] = 10.0
matplotlib.rcParams["xtick.minor.pad"] = 10.0
matplotlib.rcParams["ytick.major.pad"] = 6.0
matplotlib.rcParams["ytick.minor.pad"] = 6.0
matplotlib.rcParams["xtick.labelsize"] = 26.0
matplotlib.rcParams["ytick.labelsize"] = 26.0
matplotlib.rcParams["axes.titlesize"] = 24.0
matplotlib.rcParams["axes.labelsize"] = 28.0
matplotlib.rcParams["axes.labelpad"] = 8.0
plt.rcParams["font.size"] = 28
matplotlib.rcParams["legend.handlelength"] = 2
# matplotlib.rcParams["figure.dpi"] = 200
matplotlib.rcParams["axes.axisbelow"] = True
matplotlib.rcParams["figure.figsize"] = (13,10)
if dark:
    plt.style.use('dark_background')

gamma = 5/3.
q = 2
Geom = 4*np.pi

def fluid_eqs(xi, fluid_fields):
    a11 = lambda xi_val, D_val, V_val, P_val: V_val - (2/5)*xi_val
    a12 = lambda xi_val, D_val, V_val, P_val: D_val
    a13 = lambda xi_val, D_val, V_val, P_val: 0
    b1  = lambda xi_val, D_val, V_val, P_val: -(q/xi_val)*D_val*V_val
    a21 = lambda xi_val, D_val, V_val, P_val: 0
    a22 = lambda xi_val, D_val, V_val, P_val: V_val - (2/5)*xi_val
    a23 = lambda xi_val, D_val, V_val, P_val: 1/D_val
    b2  = lambda xi_val, D_val, V_val, P_val: (3/5)*V_val
    a31 = lambda xi_val, D_val, V_val, P_val: 0
    a32 = lambda xi_val, D_val, V_val, P_val: gamma*P_val 
    a33 = lambda xi_val, D_val, V_val, P_val: V_val - (2/5)*xi_val
    b3  = lambda xi_val, D_val, V_val, P_val: (-gamma*q/xi_val * V_val + 6/5)*P_val
    U = lambda tup: np.array([ [a11(*tup), a12(*tup), a13(*tup)],
                             [a21(*tup), a22(*tup), a23(*tup)],
                             [a31(*tup), a32(*tup), a33(*tup)],
                             ])
    B = lambda tup: np.array([[b1(*tup), b2(*tup), b3(*tup)]]).T
    U_invB = lambda tup: np.dot(np.linalg.inv(U(tup)), B(tup)).T[0]
    D, V, P = fluid_fields
    if isinstance(D, np.ndarray) and D.shape[0]>1:
        D_prime = np.zeros_like(D)
        V_prime = np.zeros_like(V)
        P_prime = np.zeros_like(P)
        for i in range(D.shape[0]):
            tup = (xi, D[i], V[i], P[i])
            D_prime[i], V_prime[i], P_prime[i] = U_invB(tup)
        return np.array([D_prime, V_prime, P_prime])
    else:
        tup = (xi, D, V, P)
        return U_invB(tup)

# RH boundary
D_1 = (gamma + 1)/(gamma-1)
V_1 = 4/(5*(gamma+1))
P_1 = 8/(25*(gamma+1))

eps = 1.0e-03 # avoid the singularity at xi=0
xi = np.linspace(1, 0+eps, 10000)

sol = solve_ivp(fluid_eqs, [1, 0+eps], [D_1, V_1, P_1], t_eval=xi,
                dense_output=True, vectorized=False,
                method="BDF")
print(sol.message)

xi = sol.t[::-1]
D, V, P = sol.y
D = D[::-1]
P = P[::-1]
V = V[::-1]


condition = np.logical_and(D>0, V>0)
condition = np.logical_and(P>0, condition)
xi = xi[condition]
D = D[condition]
V = V[condition]
P = P[condition]


k = Geom * simpson((0.5*D*V**2 + P/(gamma-1))*xi**q, x=xi)
beta = k**(-1/(3+q))

plt.title(r"$\beta = $"+f"{beta:.3f}, "+r"$\gamma = $"+f"{gamma:.2f}, "+r"$q = $"+f"{q}")
plt.plot(xi, D/D[-1], label=r"$D/D_1$")
plt.plot(xi, P/P[-1], label=r"$P/P_1$")
plt.plot(xi, V/V[-1], label=r"$V/V_1$")
plt.legend(loc="best")
plt.xlim(xmin=0+eps, xmax=1+10*eps)
plt.show()