In [ ]:
%load_ext autoreload
%autoreload 2
In [5]:
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 matplotlib.tri as tri
In [ ]:
def solver_from_2D(dom, i):
solver_func = get_local_solver_from_index_func("eruption_model", f"test_infrasound_v26_atm{dom}")
return solver_func(i)
In [18]:
# Prep interpolation grid
solver0 = solver_from_2D(1, 0)
physics = solver0.physics
base_x = np.linspace(0, 2000, 400)
base_y = np.linspace(-300, 4000, 401)
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,2,3]]
# 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]
In [19]:
def compute_ref_mapping(x_tri):
''' Compute jacobian
x_tri: triangle nodes (..., 3, 2)
Returns
Matrix collection (..., 2, 2) such that reference coordinates can be
obtained from
M @ (x - x_tri[0,:]).
'''
_a00 = x_tri[...,1:2,0:1] - x_tri[...,0:1,0:1]
_a01 = x_tri[...,2:3,0:1] - x_tri[...,0:1,0:1]
_a10 = x_tri[...,1:2,1:2] - x_tri[...,0:1,1:2]
_a11 = x_tri[...,2:3,1:2] - x_tri[...,0:1,1:2]
_M = np.zeros((*x_tri.shape[:-2], 2, 2))
_M[...,0:1,0:1] = _a11
_M[...,0:1,1:2] = -_a01
_M[...,1:2,0:1] = -_a10
_M[...,1:2,1:2] = _a00
_M /= (_a00 * _a11 - _a01 * _a10)
return _M
In [31]:
# List of file indices to read
file_index_list = np.arange(0,100,1)
In [32]:
# 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,2,3]]
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 = 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]))
In [48]:
# 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))
In [63]:
t_range = np.linspace(0, 15, T_interp.shape[0])
Compute spatial gradients numerically (using a grid)¶
In [50]:
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,...]
In [51]:
# 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
In [53]:
import matplotlib
time_idx = 10
fig, ax = plt.subplots(2,2,figsize=(8.5,9))
global_clims = (-5, 5)
levels = np.linspace(-1, 1, 100)
cf = ax[0,0].contourf(center_x, center_y, dudx[time_idx,:,:], levels=levels)
sm = plt.cm.ScalarMappable(
norm=matplotlib.colors.Normalize(vmin=cf.get_clim()[0], vmax=cf.get_clim()[1]),
cmap=matplotlib.cm.jet)
cb = plt.colorbar(sm, ax=ax[0,0])
cf = ax[0,1].contourf(center_x, center_y, dudy[time_idx,:,:], levels=levels)
sm = plt.cm.ScalarMappable(
norm=matplotlib.colors.Normalize(vmin=cf.get_clim()[0], vmax=cf.get_clim()[1]),
cmap=matplotlib.cm.jet)
cb = plt.colorbar(sm, ax=ax[0,1])
cf = ax[1,0].contourf(center_x, center_y, dvdx[time_idx,:,:], levels=levels)
sm = plt.cm.ScalarMappable(
norm=matplotlib.colors.Normalize(vmin=cf.get_clim()[0], vmax=cf.get_clim()[1]),
cmap=matplotlib.cm.jet)
cb = plt.colorbar(sm, ax=ax[1,0])
cf = ax[1,1].contourf(center_x, center_y, dvdy[time_idx,:,:], levels=levels)
sm = plt.cm.ScalarMappable(
norm=matplotlib.colors.Normalize(vmin=cf.get_clim()[0], vmax=cf.get_clim()[1]),
cmap=matplotlib.cm.jet)
cb = plt.colorbar(sm, ax=ax[1,1])
fig.tight_layout()
Lighthill analogy¶
In [54]:
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]
In [55]:
# 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
In [56]:
# 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,:]
In [57]:
# 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]
In [58]:
# 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 [61]:
fig, ax = plt.subplots(2, 4, figsize=(10.5,7))
levels = np.linspace(-.3e6, .3e6, 50) / 1e6
# t_indices = [3, 9, 36, 72]
t_indices = [0, 33, 66, 99]
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, 500)
ax[0,i].set_ylim(-125, 1000)
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 = {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(-.3e6, .3e6, 50) / 1e6
# t_indices = [3, 9, 36, 72]
t_indices = [0, 33, 66, 99]
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, 500)
ax[1,i].set_ylim(-125, 1000)
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 = {t_range[t_ind]} s")
# cb.ax.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%e'))
fig.tight_layout()
plt.draw()
plt.savefig("refblastF_lighthill_analogy_cbar.svg")
plt.savefig("refblastF_lighthill_analogy_cbar.eps")
plt.savefig("refblastF_lighthill_analogy_cbar.png", dpi=200)
