Week of 2025.05.19¶
TODO¶
- Generate cleaner derivitives of light hill stress tensor.
- Get accurate atmospheric wave as a result.
- The magnitudes of the stress tensor divergences seem sufficiently large I don't know how the accoustic pressure will not be overestimated. Once you have an animation with a slightly longer run, reach out to Eric to see if he has thoughts about what may be the issue.
- Experiment with making the atmosphere mesh smaller radius.
- Re-run the simulation with updated mesh size + increase the number of time steps to get higher frequency waves.
- Make some slides for Mario about your work
- Plot the pressure from the quail simulation at 800m and then project the continued pressure wave at 2000m.
Questions¶
- Is the very simple model the correct order of magnitude
0.0 Review eruption¶
In [3]:
%load_ext autoreload
%autoreload 2
In [19]:
from helper_code.slip_imports import *
from helper_code.helper_functions import get_local_solver_from_index_func
from helper_code.animate import animate_conduit_pressure
import helper_code.infrasound as infrasound
from matplotlib.animation import FuncAnimation
import matplotlib.pyplot as plt
plt.rcParams['animation.ffmpeg_path'] = '/opt/homebrew/bin/ffmpeg'
folder_name = "simple_eruption_model"
file_name = "tungurahua_atmosphere_1"
iterations = 1000
END_TIME = 30
C0 = 320 # m/s at 5000m
In [527]:
import matplotlib.pyplot as plt
import matplotlib.animation as animation
plt.rcParams['animation.ffmpeg_path'] = '/opt/homebrew/bin/ffmpeg'
solver_func = get_local_solver_from_index_func(folder_name, file_name)
ani = animate_conduit_pressure(solver_func, iterations=iterations, d_iterations=20, viscosity_index=1, wall_friction_index=5, max_velocity=100, max_slip=60, max_tau=0.1, max_pressure=50, max_speed_of_sound=100, max_water=20, max_density=5e3, max_fragmentation=5000, max_crystal=100, max_viscosity=1)
HTML(ani.to_html5_video())
Out[527]:
1.0 Review monopole source¶
Eric shared notes on 3D acoustics, from which I will use the relation
$$ p(r,t) = \frac{\rho \dot{Q}(t-r/c)}{4 \pi r * (2/3)} $$
See more in-depth notes in my notes from two weeks ago.
The monopole source produces unfiltered infrasound waves to the correct order of magnitude (~200 Pa) but not sufficiently precise to match the specific observed wave pattern.
In [68]:
solver_func = get_local_solver_from_index_func(folder_name, file_name)
u_vec = []
t_vec = []
R = 5 # m
for i in range(0, int(iterations), 20):
solver = solver_func(i)
momentum = solver.state_coeffs[:,:,solver.physics.get_momentum_slice()]
rho = np.sum(solver.state_coeffs[:, :, solver.physics.get_mass_slice()],axis=2,keepdims=True)
# Define velocity as momentum divided by density. "velocity" when computed as an additional state variable appears to be an absolute value.
u = momentum.ravel() / rho.ravel()
# Take only the exit velocity
u_vec.append(np.maximum(u[-1], 0))
t_vec.append(solver.time)
a_vec = np.gradient(np.array(u_vec), np.array(t_vec))
Q_dot_vec = np.pi * R**2 * a_vec # m^3/s^2
N = 100
pressure_array = np.zeros((len(t_vec), N, N))
x_low, x_upper = -500, 3000
y_low, y_upper = -1000, 2500
X = np.linspace(x_low, x_upper, N)
Y = np.linspace(y_low, y_upper, N)
# Set up the figure and axis
fig, ax = plt.subplots()
im = ax.imshow(pressure_array[0], extent=[x_low, x_upper, y_low, y_upper],
aspect='auto', cmap='jet', origin='lower', vmin=0, vmax=50)
plt.colorbar(im, label='Pressure (Pa)')
ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_title('Pressure Distribution Over Time')
# Initialization function for animation
def init():
im.set_array(np.zeros((N, N)))
return [im]
# Animation update function
def update(t_idx):
for i in range(N):
for j in range(N):
x = x_low + (x_upper - x_low) * i / N
y = y_low + (y_upper - y_low) * j / N
if infrasound.point_in_volcano(x, y):
pressure_array[int(t_idx), j, i] = 0 # Set pressure to 0 inside volcano
else:
pressure_array[int(t_idx), j, i] = infrasound.relative_pressure(t_idx, x, y, t_vec, Q_dot_vec)
im.set_array(pressure_array[int(t_idx)])
ax.set_title(f'Pressure Distribution at t={t_idx:.2f} s')
return [im]
# Create animation
ani = FuncAnimation(fig, update, init_func=init, frames=np.asarray(t_vec),
interval=200, blit=True)
display(HTML(ani.to_html5_video()))
plt.close(fig) # Clean up
In [76]:
x_idx = np.argmin(np.abs(X - 2000))
y_idx = np.argmin(np.abs(Y))
fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(10, 5), sharex=True)
axes[2].plot(t_vec, pressure_array[:, y_idx, np.argmin(np.abs(X - 2000))], label="Pressure at 2000m")
axes[1].plot(t_vec, pressure_array[:, y_idx, np.argmin(np.abs(X - 1200))], label="Pressure at 1200m")
axes[0].plot(t_vec, pressure_array[:, y_idx, np.argmin(np.abs(X - 400))], label="Pressure at 400m")
for i, ax in enumerate(axes):
ax.legend(loc='upper right')
ax.grid(True)
ax.set_ylim(-2, 4)
ax.set_xlim(0, 30)
fig.supxlabel("time(s)")
fig.supylabel("pressure (Pa)")
Out[76]:
Text(0.02, 0.5, 'pressure (Pa)')
1.1 Review Quail atmosphere model¶
- Increased atmosphere order to 1
- Reduced gride size
- Reduced radius of conduit to 5m
In [38]:
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.colors as colors
import numpy as np
import processing.mdtools as mdtools
iterations = 1000
# Assume solver2D_from is available as in the original code
solver2D_atm1 = get_local_solver_from_index_func("simple_eruption_model", "test_infrasound_v29_atm1")
# Set up figure
fig = plt.figure(figsize=(8, 6))
# Define pressure range for colormap
clims = (-1e3, 1e3)
# Set up colorbar
sm = plt.cm.ScalarMappable(
norm=colors.Normalize(vmin=clims[0], vmax=clims[1]),
cmap=plt.get_cmap()
)
cb = plt.colorbar(sm, ax=plt.gca())
cb.set_label("Pressure (Pa)")
x1, p1_0= mdtools.downsample(solver2D_atm1(0), plot_qty="Pressure")
# Animation update function
def update(frame):
plt.cla() # Clear current axis
# Fetch data for current frame
x1, p1= mdtools.downsample(solver2D_atm1(frame), plot_qty="Pressure")
# Update plot
mdtools.plot_mean(x1, p1 - p1_0, clims)
plt.title(f"Pressure Field at Sec {round(0.03 * frame)}")
return plt.gca(),
# Create animation
# Assume 100 frames for simulation; adjust based on solver data
ani = animation.FuncAnimation(
fig,
update,
frames=range(0, iterations, 50), # Adjust range based on available solver indices
interval=100, # Time between frames in milliseconds
blit=False
)
display(HTML(ani.to_html5_video()))
plt.close(fig) # Clean up
In [47]:
import matplotlib.tri as tri
from helper_code.infrasound import *
solver_2D = solver2D_atm1(0)
triangulation = tri.Triangulation(solver_2D.mesh.node_coords[...,0],
solver_2D.mesh.node_coords[...,1],
triangles=solver_2D.mesh.elem_to_node_IDs)
trifinder = triangulation.get_trifinder()
p_300 = infrasound.get_p_series(300, 0, solver2D_atm1, trifinder, iterations=1000)
p_500 = infrasound.get_p_series(500, 0, solver2D_atm1, trifinder, iterations=1000)
p_700 = infrasound.get_p_series(700, 0, solver2D_atm1, trifinder, iterations=1000)
Element ID for point (300, 0): 3082 Element ID for point (500, 0): 4468 Element ID for point (700, 0): 1153
In [59]:
fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6), sharex=True)
t = np.linspace(0, END_TIME, len(p_300))
axes[0].plot(t, p_300, label="pressure a r=300")
axes[1].plot(t, p_500, label="pressure at r=500")
axes[2].plot(t, p_700, label="pressure at r=700")
p_700_sim = np.asarray(p_300) * 300/700
t_adjusted_700 = t + 300/C0
p_2000_sim = np.asarray(p_700) * 700/2000
t_adjusted_2000 = t + 1300/C0
axes[2].plot(t_adjusted_700, p_700_sim, label="forward projected wave at r=700", c='r')
axes[3].plot(t_adjusted_2000, p_2000_sim, label="forward projected wave at r=2000", c='r')
for i, ax in enumerate(axes):
ax.legend(loc='lower left')
ax.grid(True)
ax.set_ylim(-100, 100)
ax.set_xlim(0, 30)
fig.suptitle("Pressure at various radi")
fig.supylabel("pressure (Pa)")
fig.supxlabel("time (s)")
print(f"Max pressure at 2000m is {max(p_2000_sim)} Pa")
Max pressure at 2000m is 43.81281031234539 Pa
2.0 Divergence and Double divergence of lighthill stress tensor¶
In [450]:
import helper_code.lighthill as lighthill
import numpy as np
import matplotlib.tri as tri
def solver_from_2D(dom, i):
solver_func = get_local_solver_from_index_func("simple_eruption_model", f"test_infrasound_v29_atm{dom}")
return solver_func(i)
# Prep interpolation grid
solver0 = solver_from_2D(1, 0)
physics = solver0.physics
base_x = np.linspace(0, 300, 200)
base_y = np.linspace(-200, 800, 200)
mg_x, mg_y = np.meshgrid(base_x, base_y, indexing="xy")
solver0_list = [solver_from_2D(dom_idx, 0) for dom_idx in [1]]
# Compute workload partition
ind_partition = [np.where(tri.Triangulation(
solver.mesh.node_coords[...,0],
solver.mesh.node_coords[...,1],
triangles=solver.mesh.elem_to_node_IDs).get_trifinder()(mg_x.ravel(), mg_y.ravel()) != -1)[0]
for solver in solver0_list]
# List of file indices to read
file_index_list = np.arange(0,iterations,5)
In [451]:
# Allocate union (joining all times) U, in spatially-flattened shape
U_union = np.nan * np.empty((file_index_list.size, *mg_x.ravel().shape, 8+physics.NDIMS))
for time_idx, file_idx in enumerate(file_index_list):
# Load solvers for given time_idx
solver_list = [solver_from_2D(dom_idx, file_idx) for dom_idx in [1]]
for solver, _index_partition in zip(solver_list, ind_partition):
# Identify physical position (x, y) of points to interpolate at with shape (npoints, 2)
_phys_pos = np.stack(
(mg_x.ravel()[_index_partition],
mg_y.ravel()[_index_partition]),
axis=1)
# Identify element indices for all points to interpolate at
elt_indices = tri.Triangulation(
solver.mesh.node_coords[...,0],
solver.mesh.node_coords[...,1],
triangles=solver.mesh.elem_to_node_IDs).get_trifinder()(_phys_pos[:,0], _phys_pos[:,1])
# Identify element node coordinates
x_tri = solver.mesh.node_coords[solver.mesh.elem_to_node_IDs]
# Compute global physical-to-reference coordinate mapping
ref_mapping = lighthill.compute_ref_mapping(x_tri)
# Interpolate for each point using the correct element, writing to correct index in global U array
for (write_idx, x_point, ie) in zip(_index_partition, _phys_pos, elt_indices):
# For element containing point, compute reference coordinate of sampling point
ref_coords_loc = np.einsum("ij, j -> i",
ref_mapping[ie,...],
x_point - x_tri[ie,0,:])
# Evaluate basis at reference coordinate
#U_union[time_idx,write_idx,:] = (solver.state_coeffs[ie,0,:] * (1 - ref_coords_loc[0] - ref_coords_loc[1]))
U_union[time_idx,write_idx,:] = (solver.state_coeffs[ie,0,:] * (1 - ref_coords_loc[0] - ref_coords_loc[1]) + solver.state_coeffs[ie,1,:] * ref_coords_loc[0] + solver.state_coeffs[ie,2,:] * ref_coords_loc[1])
In [452]:
# Evaluate temperature using interpolated state, migrate to meshgrid shape (time_indices, mg_x.shape[0], mg_x.shape[1])
T_interp = np.reshape(physics.compute_variable("Temperature", U_union).squeeze(axis=2), (file_index_list.size, *mg_x.shape))
rho = U_union[...,0:3].sum(axis=-1, keepdims=True)
yM = np.reshape(U_union[...,2:3] / rho, (file_index_list.size, *mg_x.shape))
t_range = np.linspace(0, 30, T_interp.shape[0])
In [453]:
rho = U_union[...,0:3].sum(axis=-1, keepdims=True)
u = U_union[:,:,3:4] / rho
v = U_union[:,:,4:5] / rho
# Pull rho, u, v
mg_u = np.reshape(u, (file_index_list.size, *mg_x.shape))
mg_v = np.reshape(v, (file_index_list.size, *mg_x.shape))
mg_T = np.reshape(physics.compute_variable("Temperature", U_union).squeeze(axis=2), (file_index_list.size, *mg_x.shape))
mg_p = np.reshape(physics.compute_variable("Pressure", U_union).squeeze(axis=2), (file_index_list.size, *mg_x.shape))
mg_rho = np.reshape(rho, (file_index_list.size, *mg_x.shape))
mg_c = np.reshape(physics.compute_variable("SoundSpeed", U_union).squeeze(axis=2), (file_index_list.size, *mg_x.shape))
mg_p0 = mg_p[0,...]
/Users/paxton/git/quail_volcano/src/physics/multiphasevpT/multiphasevpT.py:279: RuntimeWarning: invalid value encountered in sqrt varq = np.sqrt(get_Gamma()*get_pressure()/(arhoA+arhoWv+arhoM)*get_Psi1())
In [454]:
# Grid dimensions
dx = np.diff(mg_x[0:1,:], axis=1)
dy = np.diff(mg_y[:,0:1], axis=0)
# Grid-center coordinates
center_x = 0.5 * (mg_x[:,1:] + mg_x[:,:-1])[:-1,:]
center_y = 0.5 * (mg_y[1:,:] + mg_y[:-1,:])[:,:-1]
# Grid-center differentiation
u_foldx = (0.5 * (mg_u[:,:,1:] + mg_u[:,:,:-1]))
u_foldy = (0.5 * (mg_u[:,1:,:] + mg_u[:,:-1,:]))
dudy = np.diff(u_foldx, axis=1) / dy
dudx = np.diff(u_foldy, axis=2) / dx
v_foldx = (0.5 * (mg_v[:,:,1:] + mg_v[:,:,:-1]))
v_foldy = (0.5 * (mg_v[:,1:,:] + mg_v[:,:-1,:]))
dvdy = np.diff(v_foldx, axis=1) / dy
dvdx = np.diff(v_foldy, axis=2) / dx
mg_c0 = mg_c[0,...]
# Grid-center coordinates
center_x = 0.5 * (mg_x[:,1:] + mg_x[:,:-1])[:-1,:]
center_y = 0.5 * (mg_y[1:,:] + mg_y[:-1,:])[:,:-1]
# Interior grid-center coordinates
int_x = center_x[1:-1,1:-1]
int_y = center_y[1:-1,1:-1]
Double divergence nonlinear¶
In [455]:
# Evaluate matrix T
T00_nonlinear = mg_p - mg_c0*mg_c0 * mg_rho
T01_nonlinear = np.zeros_like(mg_rho)
T10_nonlinear = np.zeros_like(mg_rho)
T11_nonlinear = mg_p - mg_c0*mg_c0 * mg_rho
# Mixed derivative (symmetric)
d2T01 = (np.diff(np.diff(T01_nonlinear, axis=-1), axis=-2) / dx) / dy
# d^2(T00)/dx^2
T00_foldy = 0.5 * (T00_nonlinear[:,1:,:] + T00_nonlinear[:,:-1,:])
d2T00 = T00_foldy[:,:,3:] - T00_foldy[:,:,2:-1] - T00_foldy[:,:,1:-2] + T00_foldy[:,:,0:-3]
# d^2(T11)/dy^2
T11_foldx = 0.5 * (T11_nonlinear[:,:,1:] + T11_nonlinear[:,:,:-1])
d2T11 = T11_foldx[:,3:,:] - T11_foldx[:,2:-1,:] - T11_foldx[:,1:-2,:] + T11_foldx[:,0:-3,:]
# Evaluate Lighthill analogy source at cell center, for interior cells (discard boundary cells)
dijTij_nonlinear = d2T00[:,1:-1,:] + 2 * d2T01[:,1:-1,1:-1] + d2T11[:,:,1:-1]
Double divergence inertial¶
In [456]:
# Evaluate matrix T
T00_inertial = mg_rho * mg_u * mg_u
T01_inertial = mg_rho * mg_u * mg_v
T10_inertial = mg_rho * mg_v * mg_u
T11_inertial = mg_rho * mg_v * mg_v
# Mixed derivative (symmetric)
d2T01 = (np.diff(np.diff(T01_inertial, axis=-1), axis=-2) / dx) / dy
# d^2(T00)/dx^2
T00_foldy = 0.5 * (T00_inertial[:,1:,:] + T00_inertial[:,:-1,:])
d2T00 = T00_foldy[:,:,3:] - T00_foldy[:,:,2:-1] - T00_foldy[:,:,1:-2] + T00_foldy[:,:,0:-3]
# d^2(T11)/dy^2
T11_foldx = 0.5 * (T11_inertial[:,:,1:] + T11_inertial[:,:,:-1])
d2T11 = T11_foldx[:,3:,:] - T11_foldx[:,2:-1,:] - T11_foldx[:,1:-2,:] + T11_foldx[:,0:-3,:]
# Evaluate Lighthill analogy source at cell center, for interior cells (discard boundary cells)
dijTij_inertial = d2T00[:,1:-1,:] + 2 * d2T01[:,1:-1,1:-1] + d2T11[:,:,1:-1]
In [461]:
import matplotlib
fig, ax = plt.subplots(2, 4, figsize=(10.5,7))
levels = np.linspace(-6e4, 6e4, 50) / 1e6
# t_indices = [3, 9, 36, 72]
t_indices = [0, 40, 80, 120]
for i, t_ind in enumerate(t_indices):
ax[0,i].set_facecolor("black")
cf = ax[0,i].contourf(int_x, int_y, np.clip(dijTij_inertial[t_ind,:,:] / 1e6, levels[0], levels[-1]), levels=levels, cmap=matplotlib.cm.PiYG)
cb = fig.colorbar(cf, label=r"$\partial_i \partial_j (\rho u_i u_j)$ (MPa / m${}^2$)",
# format=matplotlib.ticker.FormatStrFormatter('%.3e')
)
ax[0,i].set_xlim(0, 200)
ax[0,i].set_ylim(-100, 500)
#ax[0,i].set_aspect('equal')
ax[0,i].set_xlabel("$r$ (m)")
ax[0,i].set_ylabel("$z$ (m)")
curr_ax = ax[0,i]
for item in ([curr_ax.xaxis.label, curr_ax.yaxis.label] +
curr_ax.get_xticklabels() + curr_ax.get_yticklabels()):
item.set_fontsize(11)
cb.set_ticks(np.linspace(levels[0], levels[-1], 7))
ax[0,i].set_title(f"t = {round(t_range[t_ind])} s")
# cb.ax.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%e'))
# fig.tight_layout()
# fig, ax = plt.subplots(1, 3, figsize=(8.5,3.5))
levels = np.linspace(-1e5, 1e5, 50) / 1e6
# t_indices = [3, 9, 36, 72]
t_indices = [0, 40, 80, 120]
for i, t_ind in enumerate(t_indices):
ax[1,i].set_facecolor("black")
cf = ax[1,i].contourf(int_x, int_y, np.clip(dijTij_nonlinear[t_ind,:,:] / 1e6, levels[0], levels[-1]), levels=levels, cmap=matplotlib.cm.PiYG)
cb = fig.colorbar(cf, label=r"$\nabla^2 (p - c_0^2 \rho)$ (MPa / m${}^2$)",
# format=matplotlib.ticker.FormatStrFormatter('%.3e')
)
ax[1,i].set_xlim(0, 200)
ax[1,i].set_ylim(-100, 500)
#ax[1,i].set_aspect('equal')
ax[1,i].set_xlabel("$r$ (m)")
ax[1,i].set_ylabel("$z$ (m)")
curr_ax = ax[1,i]
for item in ([curr_ax.xaxis.label, curr_ax.yaxis.label] +
curr_ax.get_xticklabels() + curr_ax.get_yticklabels()):
item.set_fontsize(11)
cb.set_ticks(np.linspace(levels[0], levels[-1], 7))
ax[1,i].set_title(f"t = {round(t_range[t_ind])} s")
# cb.ax.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%e'))
fig.tight_layout()
plt.draw()
plt.savefig(f'{BASE_PATH}/volcano_sims/notebooks/charts/lighthill_tensor_double_divergence.png')
Single divergence nonlinear¶
In [458]:
Nt, Nx, Ny = mg_p.shape
# Central differences using np.diff (second difference)
d_p_dx = np.diff(mg_p, n=2, axis=1) / (2 * dx[0, 0])
d_p_dy = np.diff(mg_p, n=2, axis=2) / (2 * dy[0, 0])
d_rho_dx = np.diff(mg_rho, n=2, axis=1) / (2 * dx[0, 0])
d_rho_dy = np.diff(mg_rho, n=2, axis=2) / (2 * dy[0, 0])
# Pad arrays to restore original shape (Nx, Ny)
diTij_nonlinear = np.zeros((2, Nt, Nx, Ny))
d_c0_dx0 = np.diff(mg_c0, n=2, axis=0) / (2 * dx[0, 0])
d_c0_dx1 = np.diff(mg_c0, n=2, axis=1) / (2 * dy[0, 0])
diTij_nonlinear[0, :, 1:-1, :] = d_p_dx - mg_c0[1:-1, :]**2 * d_rho_dx - 2 * mg_c0[1:-1, :] * mg_rho[:, 1:-1, :] * d_c0_dx0
diTij_nonlinear[1, :, :, 1:-1] = d_p_dy - mg_c0[:, 1:-1]**2 * d_rho_dy - 2 * mg_c0[:, 1:-1] * mg_rho[:, :, 1:-1] * d_c0_dx1
In [459]:
clims = [-1e4, 1e4]
fig, ax = plt.subplots(1, 2)
img = ax[0].imshow(diTij_nonlinear[0,0], aspect='auto', vmin=clims[0], vmax=clims[1], extent=[int_x[0,0], int_x[0,-1], int_y[-1,0], int_y[0, 0]], cmap=matplotlib.cm.PiYG)
img2 = ax[1].imshow(diTij_nonlinear[1,0], aspect='auto', vmin=clims[0], vmax=clims[1], extent=[int_x[0,0], int_x[0,-1], int_y[-1,0], int_y[0, 0]], cmap=matplotlib.cm.PiYG)
plt.colorbar(img)
plt.colorbar(img2)
# Adjust spacing between subplots
plt.subplots_adjust(hspace=0.6)
plt.suptitle("Divergence of Lighthill Tensor")
ax[0].set_title('r-component')
ax[0].set_ylabel('height (m)')
ax[0].invert_yaxis()
ax[1].set_title('y-component')
ax[1].set_xlabel('radius (m)')
ax[1].invert_yaxis()
# Update function for animation
def update(i):
img.set_array(diTij_nonlinear[0,i])
img2.set_array(diTij_nonlinear[1,i])
return [img, img2]
# Create animation
ani = FuncAnimation(fig, update, frames=np.arange(len(t_range)),
interval=50, blit=True)
display(HTML(ani.to_html5_video()))
ani.save(f'{BASE_PATH}/volcano_sims/notebooks/animations/lighthill_tensor_divergence.mp4', writer='ffmpeg', fps=30, bitrate=1800)
plt.close(fig) # Clean up
3.0 Volume integral approach¶
I am calculating the volume integral over the half-cylindrical region bounded by $r < 400$ and $y > -100$ and $y<500$.
As we can see in the above graphics, the value of the double divergence of the lighthill stress tensor falls off outside those bounds.
The equation we solve below is:
$$ p'(x,t) = \frac{1}{4\pi} \int_v \frac{1}{|x - y|} \frac{\partial^2 T_{ij}}{\partial y_i \partial y_j} (y, t - \frac{|x-y|}{c_0}) dy^3 $$
To solve for pressure as a time series, we can just iterate over the integral for all time steps. We need only calculate the value of $\frac{\partial^2 T_{ij}}{\partial y_i \partial y_j}$ when $t - \frac{|x-y|}{c_0}$ is greater than 0 and less than the end of the simulation.
- The $4 \pi$ in the denominator is technically incorrect for our geometry.
- Once I feel confident about the soundness of my code, I should decrease the mesh size and time sampling size to get a more precise result.
My implementation can be found here.
In [460]:
import helper_code.lighthill as lighthill
source = dijTij_inertial + dijTij_nonlinear
print(source.shape)
x_obs = (2000, -1000, 0)
c0 = 320 # m/s roughly the speed of sound in air at 5000m in elevation
atmosphere_solid_angle = 4 * np.pi * (2/3)
X = int_x[0]
Y = int_y[:,0]
X = np.linspace(-200, 200, 20)
Y = np.linspace(0, 500, 20)
Z = np.linspace(-200, 200, 20)
points = (t_range, int_y[:,0], int_x[0])
pressure_obs = lighthill.calculate_pressure_as_volume_integral(X, Y, Z, file_index_list, x_obs, points, source, t_range, int_x[0][-1], c0)
(200, 197, 197) DV size is 11663.507799970843 Finished time index 0 of 200 Finished time index 10 of 200 Finished time index 20 of 200 Finished time index 30 of 200 Finished time index 40 of 200 Finished time index 50 of 200 Finished time index 60 of 200 Finished time index 70 of 200 Finished time index 80 of 200 Finished time index 90 of 200 Finished time index 100 of 200 Finished time index 110 of 200 Finished time index 120 of 200 Finished time index 130 of 200 Finished time index 140 of 200 Finished time index 150 of 200 Finished time index 160 of 200 Finished time index 170 of 200 Finished time index 180 of 200 Finished time index 190 of 200 Number of contributions: 1203028 r average: 2363.7261661976113 Max source value: [24689331.80069273]
In [469]:
plt.plot(t_range, lighthill.highpass(pressure_obs), label="simulated data based on volume integral")
plt.plot(tr0_filtered.times(), tr0_filtered.data, label="observed data")
plt.xlabel('Time (s)')
plt.ylabel('Acoustic Pressure (Pa)')
plt.title('Filtered Acoustic Pressure at Observer based on Lighthill analogy')
plt.legend()
Out[469]:
<matplotlib.legend.Legend at 0x11ebc7250>
3.1 validate model with known source¶
$$ G(\bf{x}, \bf{y}, t) = \frac{\delta(t - \frac{|\bf{x} - \bf{y}|}{c_0})}{4 \pi | \bf{x} - \bf{y}|} $$
So the pressure at any given point can be expressed as
$$ \begin{align} p(\bf{x}, t) =& \int \int_V G(\bf{x}, \bf{y}, t - \tau) S(\bf{y}, \tau) d^3y d \tau \\ \end{align} $$
In this case, because the source term is uniform across the volume $R < 100$ and $Y>0$ and $Y< 100$.
To find the max pressure we would expect, we write the equation:
$$ p_{max} = \frac{1}{4 \pi} \int_V \frac{1}{|x_{obs} - y|} dy $$
That can roughly be approximated as
$$ \begin{align} p_{max} =& \frac{1}{4 \pi} \frac{V}{r_{avg}} \\ p_{max} =& \frac{1}{4 \pi} \frac{\pi * 100^2 * 200}{2300} \\ p_{max} =& 220 [Pa] \end{align} $$
In the below example, I show that numerically I arrive at a solution with a max pressure of about 190.
In [442]:
n = len(t_range) # Length of the vector
mu = n // 4 # Center of the pulse is 1/4 through the time series
sigma = 5 # Standard deviation (controls pulse width)
nx = len(int_x[0])
ny = len(int_y[:,0])
print(f"nx {nx} and ny {ny}")
# Create index array
x = np.arange(n)
# Compute Gaussian pulse
pulse = np.exp(-((x - mu) ** 2) / (2 * sigma ** 2))
source_gaussian = np.zeros(dijTij_inertial.shape)
nan_mask = np.isnan(dijTij_inertial)
for i in range(nx):
for j in range(ny):
if int_y[j, 0] > 0 and int_y[j,0] < 200 and int_x[0,i] < 100:
source_gaussian[:, j, i] = pulse
source_gaussian[nan_mask] = np.nan
clims = [-1, 1]
fig, ax = plt.subplots()
img = ax.imshow(source_gaussian[0], cmap=matplotlib.cm.PiYG, aspect='auto', vmin=clims[0], vmax=clims[1], extent=[int_x[0,0], int_x[0,-1], int_y[-1,0], int_y[0, 0]])
plt.gca().invert_yaxis()
plt.colorbar(img)
ax.set_title('Uniform Gaussian Pulse')
ax.set_xlabel('radius (m)')
ax.set_ylabel('height (m)')
# Update function for animation
def update(i):
img.set_array(source_gaussian[i])
return [img]
# Create animation
ani = FuncAnimation(fig, update, frames=np.arange(len(t_range)),
interval=50, blit=True)
display(HTML(ani.to_html5_video()))
plt.close(fig) # Clean up
nx 197 and ny 197
In [ ]:
X = np.linspace(-200, 100, 20)
Y = np.linspace(0, 200, 20)
Z = np.linspace(-100, 100, 20)
points = (t_range, int_y[:,0], int_x[0])
pressure_test = lighthill.calculate_pressure_as_volume_integral(X, Y, Z, file_index_list, x_obs, points, source_gaussian, t_range, int_x[0][-1], c0)
DV size is 1749.5261699956259 Finished time index 0 of 100 Finished time index 10 of 100 Finished time index 20 of 100 Finished time index 30 of 100 Finished time index 40 of 100 Finished time index 50 of 100 Finished time index 60 of 100 Finished time index 70 of 100 Finished time index 80 of 100 Finished time index 90 of 100 Number of contributions: 506840 r average: 2326.6972476922006 Max source value: [0.999992234645567]
In [ ]:
r_avg = 2300 #m
retarded_time = t_range - r_avg/c0
# Plot acoustic pressure
plt.plot(t_range, pressure_test, label="numerical result plotted against retarded time")
plt.plot(t_range, np.ones(t_range.shape)*220, label="Analytically predicted max pressure (220 Pa)")
plt.xlabel('Time (s)')
plt.ylabel('Acoustic Pressure (Pa)')
plt.title('Acoustic Pressure at Observer')
plt.grid(True)
plt.legend()
plt.show()
print(f"Max pressure {max(pressure_test)}")
Max pressure 211.2469422199135
4.0 Surface integral flux approach¶
Eric suggested I solved this problem as a surface integral. Doing so would allow us to remove one spatial dimension at the expense of slightly more complicated math.
Previously, we had expressed the pressure as the following integral.
$$ \begin{align} p'(x,t) = \frac{1}{4\pi} \int_v \frac{1}{|x - y|} \frac{\partial}{\partial y_i} \frac{\partial T_{ij}}{\partial y_j} (y, t - \frac{|x-y|}{c_0}) dy^3 \end{align} $$
Applying the divergence theorm, we are able to rewrite this volume integral as a surface intgral.
$$ \begin{align} p'(x,t) = \frac{1}{4\pi} \int_S \frac{1}{|x - y|} \frac{\partial T_{ij}}{\partial y_j} (y, t - \frac{|x-y|}{c_0}) \cdot \hat{n} dy^2 \end{align} $$
Let's integrate over the sphere with a radius $a=100m$. Let's review a couple aspects of spherical coordinates.
$$ \begin{align} dS =& a^2 \sin \phi d \phi d \theta \\ x =& a \sin \phi \cos \theta \\ y =& a \sin \phi \sin \theta \\ z =& a \cos \phi \end{align} $$
So we should be able to rewrite the surface integral as follows:
$$ p'(x, t) = \frac{1}{4 \pi} \int_0^{\pi/3} \int_0^{2\pi} \frac{1}{|x - y|} \frac{\partial T_{ij}}{\partial y_j} n_j a^2 \sin \phi d \phi d \theta $$
In [471]:
points = (range(2), t_range, mg_y[:,0], mg_x[0])
pressure_from_surface = np.zeros(len(file_index_list))
for i, time_idx in enumerate(file_index_list):
pressure_from_surface[i] = lighthill.calculate_surface_integral(diTij_nonlinear, 50, i, t_range, x_obs, c0, points)
if i % 10 == 0:
print(f"Completed step of index {i}")
Completed step of index 0 Completed step of index 10 Completed step of index 20 Completed step of index 30 Completed step of index 40
/var/folders/f7/g7y34v812n5_58yvtkt3mrw40000gn/T/ipykernel_97421/1753142958.py:6: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) pressure_from_surface[i] = lighthill.calculate_surface_integral(diTij_nonlinear, 50, i, t_range, x_obs, c0, points)
Completed step of index 50 Completed step of index 60 Completed step of index 70 Completed step of index 80 Completed step of index 90 Completed step of index 100 Completed step of index 110 Completed step of index 120 Completed step of index 130 Completed step of index 140 Completed step of index 150 Completed step of index 160 Completed step of index 170 Completed step of index 180 Completed step of index 190
In [524]:
#plt.plot(t_range, pressure)
plt.plot(t_range, lighthill.highpass(pressure_from_surface))
plt.ylabel("Pressure at observer (Pa)")
plt.xlabel("time (s)")
Out[524]:
Text(0.5, 0, 'time (s)')
A1. Observed infrasound data¶
In [ ]:
import matplotlib.pyplot as plt
import numpy as np
import obspy
from obspy import read, Stream, UTCDateTime
from obspy.io.mseed.util import get_record_information
file_path = "/Users/paxton/git/volcano_sims/infrasound_data/8A.HIGH.2013-07-14.mseed"
st = read(file_path)
tr0 = st[0]
tr0_filtered = tr0.copy()
freqmin = 1
freqmax = 20
tr0_filtered.filter("bandpass", freqmin=freqmin, freqmax=freqmax, corners=4, zerophase=True)
tr0_filtered.data = tr0_filtered.data * 1.589e-6 / 23e-6
t1 = UTCDateTime("2013-07-14T11:46:30")
t2 = UTCDateTime("2013-07-14T11:46:45")
tr0_filtered.trim(t1, t2)
tr0_filtered.plot()
Out[ ]: