Skip to content

Instantly share code, notes, and snippets.

@dutta-alankar

dutta-alankar/sedov-taylor.html

Created February 17, 2025 16:12
Show Gist options
  • Save dutta-alankar/70145306cbdece19103dd511fa9b94c9 to your computer and use it in GitHub Desktop.
Save dutta-alankar/70145306cbdece19103dd511fa9b94c9 to your computer and use it in GitHub Desktop.
Download ZIP
Sedov-Taylor blastwave solution
Raw
sedov-taylor.html
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Sedov-Taylor blastwave solver</title>
<script src="https://cdn.jsdelivr.net/pyodide/v0.25.0/full/pyodide.js"></script>
<style>
:root {
--bg-color: white;
--text-color: black;
--btn-bg: #ddd;
--btn-text: black;
}
[data-theme="dark"] {
--bg-color: #1e1e1e;
--text-color: white;
--btn-bg: #444;
--btn-text: white;
}
body {
background-color: var(--bg-color);
color: var(--text-color);
font-family: Arial, sans-serif;
text-align: center;
transition: background 0.3s ease-in-out, color 0.3s ease-in-out;
}
.container {
display: flex;
flex-direction: column;
align-items: center;
gap: 15px;
max-width: 600px;
margin: auto;
}
.controls {
display: flex;
flex-direction: column;
align-items: center;
gap: 10px;
width: 100%;
}
button {
background-color: var(--btn-bg);
color: var(--btn-text);
border: none;
padding: 10px 15px;
margin: 10px;
cursor: pointer;
transition: background 0.3s, transform 0.2s;
}
button:hover {
transform: scale(1.05);
}
img {
border: 0px solid var(--text-color);
opacity: 0;
transition: opacity 0.5s ease-in-out;
}
.loading {
font-size: 1.5em;
font-weight: bold;
animation: blink 1s infinite;
}
@keyframes blink {
0% { opacity: 1; }
50% { opacity: 0.5; }
100% { opacity: 1; }
}
</style>
</head>
<body>
<h1>Sedov-Taylor blastwave solver</h1>
<p id="pyodide-loading" class="loading">Loading Pyodide...</p>
<div class="container" id="main-container" style="display:none;">
<button onclick="toggleTheme()">Toggle Light/Dark mode</button>
<div class="controls">
<label for="gamma-slider">Adiabatic index Gamma: <span id="gamma-value">1.66667</span></label>
<input type="range" id="gamma-slider" min="0.1" max="2.0" step="0.1" value="1.6667" disabled>
<label>Geometry</label>
<select id="geom-select" disabled>
<option value="2">Spherical</option>
<option value="1">Cylindrical</option>
<option value="0">Cartesian</option>
</select>
</div>
<button id="solve-btn" onclick="solveODEFirstTime()" disabled>Solve</button>
<p id="loading" class="loading" style="display:none;">Solving ODE...</p>
<img id="plot" src="" alt="ODE solution will appear here">
<button id="export-btn" onclick="exportImage()" style="display:none;">Export Image</button>
</div>
<script>
let pyodideReady = false;
let firstSolve = false;
function detectSystemTheme() {
return window.matchMedia("(prefers-color-scheme: dark)").matches ? "dark" : "light";
}
function applyTheme(theme) {
document.documentElement.setAttribute("data-theme", theme);
localStorage.setItem("theme", theme);
}
function toggleTheme() {
let currentTheme = document.documentElement.getAttribute("data-theme");
let newTheme = currentTheme === "dark" ? "light" : "dark";
applyTheme(newTheme);
solveODE();
}
function loadTheme() {
let savedTheme = localStorage.getItem("theme") || detectSystemTheme();
applyTheme(savedTheme);
}
async function loadPyodideAndPackages() {
console.log("Loading Pyodide...");
window.pyodide = await loadPyodide();
await pyodide.loadPackage(["numpy", "scipy", "matplotlib"]);
console.log("Pyodide and required packages installed.");
pyodideReady = true;
document.getElementById("pyodide-loading").style.display = "none";
document.getElementById("main-container").style.display = "flex";
document.getElementById("gamma-slider").disabled = false;
document.getElementById("geom-select").disabled = false;
document.getElementById("solve-btn").disabled = false;
}
async function solveODEFirstTime() {
document.getElementById("solve-btn").style.display = "none";
firstSolve = true;
solveODE();
}
async function solveODE() {
if (!pyodideReady) {
console.log("Pyodide not yet loaded. Waiting...");
return;
}
if (!firstSolve) {
console.log("Hit the solve button for the first time ...");
let gamma = parseFloat(document.getElementById("gamma-slider").value);
let q = parseFloat(document.getElementById("geom-select").value);
document.getElementById("gamma-value").innerText = gamma;
let theme = document.documentElement.getAttribute("data-theme");
console.log("q = ", q);
console.log("gamma = ", gamma);
return;
}
let gamma = parseFloat(document.getElementById("gamma-slider").value);
let q = parseFloat(document.getElementById("geom-select").value);
document.getElementById("gamma-value").innerText = gamma;
let theme = document.documentElement.getAttribute("data-theme");
console.log("q = ", q);
console.log("gamma = ", gamma);
document.getElementById("loading").style.display = "block";
document.getElementById("plot").style.opacity = "0";
let pythonCode = `
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp, simpson
import io, base64
dark = "${theme}" == "dark"
plt.style.use("dark_background" if dark else "default")
## 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["axes.axisbelow"] = True
matplotlib.rcParams["figure.figsize"] = (13,10)
gamma = np.abs(${gamma})
q = np.abs(${q})
if q==0:
Geom =0
elif q==1:
Geom = 2*np.pi
else:
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.figure(figsize=(13,10))
plt.title(r"$\\beta = $"+f"{beta:.3f}, "+r"$\\gamma = $"+f"{gamma:.2f}, "+r"$q = $"+f"{q}")
plt.text(0.45, 0.9, r"$R_s(t) = \\beta \\left(\\frac{E_0 t^2}{\\rho_0}\\right)^{\\frac{1}{3+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)
buf = io.BytesIO()
plt.savefig(buf, format="png", transparent=True)
buf.seek(0)
"data:image/png;base64," + base64.b64encode(buf.read()).decode()
`;
let result = await pyodide.runPythonAsync(pythonCode);
document.getElementById("loading").style.display = "none";
document.getElementById("plot").src = result;
document.getElementById("plot").style.opacity = "1";
document.getElementById("export-btn").style.display = "block";
}
function exportImage() {
let img = document.getElementById("plot");
let link = document.createElement("a");
link.href = img.src;
link.download = "ODE_Solution.png";
document.body.appendChild(link);
link.click();
document.body.removeChild(link);
}
document.getElementById("gamma-slider").addEventListener("input", solveODE);
document.getElementById("geom-select").addEventListener("change", solveODE);
loadTheme();
loadPyodideAndPackages();
</script>
</body>
</html>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment