In [8]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
anim = animate2d(Q, x, y, T, fps = 24, frames = 144)
display(anim)
import cmocean as cmo
import time, sys
import numpy as np
import scipy.constants as spc
from scipy import signal as sg
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation
from IPython.display import display, Math, HTML
from numpy.linalg import inv
from matplotlib.animation import FuncAnimation
from scipy.ndimage import convolve, generate_binary_structure
plt.rcParams.update(
{
'mathtext.fontset':'cm',
'font.family':'serif',
'font.sans-serif':'Times New Roman',
'figure.figsize':[9,4],
'figure.titlesize':22,
"figure.dpi":90,
'savefig.dpi':90,
'axes.titlesize':20,
'axes.labelsize':18,
'axes.titley': 1.0,
'axes.titlepad': 5.0,
'axes.edgecolor':'black',
'axes.grid': False,
'grid.alpha': .5,
'xtick.labelsize':14,'ytick.labelsize':14,
'xtick.major.size':6,'ytick.major.size':6,
'xtick.major.width':1.25, 'ytick.major.width':1.25,
'xtick.direction':'inout','ytick.direction':'inout',
'xtick.top':False, 'ytick.right':False,
'legend.title_fontsize':14, 'legend.fontsize':14,
'legend.borderaxespad': 1, 'legend.borderpad': 0.5,
'legend.framealpha': 1,
'legend.handleheight': 0.5, 'legend.handlelength': 2.0, 'legend.handletextpad':0.5,
'legend.labelspacing': 0.25,
'legend.fancybox':False,
'legend.edgecolor': '0',
'legend.frameon': True,
'legend.markerscale': 1.25,
'animation.embed_limit':2**128,
'animation.html': 'jshtml'
}
)
pi = spc.pi
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
default_plots()
def animate2d(u, x, y, t, fps = 24, frames = 144):
X, Y = np.meshgrid(x, y)
t_steps = int(len(t)-1)
i_frames = np.arange(0, t_steps, int(t_steps/frames))
fig = plt.figure(figsize = (6,6), constrained_layout = True)
ax = fig.add_subplot(111, projection = '3d')
surf = ax.plot_surface(X, Y, u[0], alpha = 1)
def animate(i):
ax.clear()
surf = ax.plot_surface(X, Y, u[i], cmap = cmo.cm.dense, alpha = 1)
ax.axis(False)
ax.set_xlim(X.min(), X.max())
ax.set_ylim(Y.min(), Y.max())
ax.set_zlim(u[0].min(), u[0].max())
return surf,
anim = FuncAnimation(fig, animate, interval = int(1e3/fps), frames = i_frames, blit = False)
plt.close()
return anim
def wave2D(q0, p0, x, y, mask, c = 1, t0 = 0, tf = 20, g = 1, boundary = 'periodic'):
def fix(u):
u[0,:] = u[:,0] = u[-1,:] = u[:,-1] = 0
return
dx = x[1] - x[0]
dy = y[1] - y[0]
dt = (g/c)*(dx+dy)/2
δx = (c/dx)**2
δy = (c/dy)**2
steps = int((tf-t0)/dt)
t = np.linspace(t0, tf, steps+1)
q = np.zeros((steps+1,) + q0.shape)
p = np.zeros((steps+1,) + p0.shape)
q[0] = q0
fix(q.T)
p[0] = p0
if boundary == 'periodic':
for i in range(steps):
q_R, q_L = np.roll(q[i],(0,-1),(0,1)), np.roll(q[i],(0,1),(0,1))
q_T, q_B = np.roll(q[i],(1,0),(0,1)), np.roll(q[i],(-1,0),(0,1))
Dx = (q_R + q_L - 2*q[i]) * δx
Dy = (q_T + q_B - 2*q[i]) * δy
p[i+1] = p[i] + dt * (Dx + Dy)
q[i+1] = q[i] + dt * p[i+1]
elif boundary == 'fixed':
Δ1 = np.array([[0, 0, 0], [1, -2, 1], [0, 0, 0]]) * δx
Δ2 = np.array([[0, 1, 0], [0, -2, 0], [0, 1, 0]]) * δy
for i in range(steps):
fix(q.T)
q[i] *= mask
dp_dt1 = convolve(q[i], Δ1, mode = 'constant')[1:-1,1:-1]
dp_dt2 = convolve(q[i], Δ2, mode = 'constant')[1:-1,1:-1]
p[i+1,1:-1,1:-1] = p[i,1:-1,1:-1] + dt * (dp_dt1 + dp_dt2)
q[i+1,1:-1,1:-1] = q[i,1:-1,1:-1] + dt * p[i+1,1:-1,1:-1]
fix(q.T)
return q, p, t
Lx = 6
Nx = 128
dx = Lx/(Nx-1)
Ly = Lx
Ny = Nx
dy = dx
x = np.arange(-1, Nx+1)*dx
y = np.arange(-1, Ny+1)*dy
a = 3/2
X, Y = np.meshgrid(x,y)
R = 1
c_mask = ((np.sqrt((X-Lx/2)**2+(Y-Ly/2)**2) <= 3*R)).astype(float)
q0 = np.sqrt(4/(Lx*Ly))*np.sin(3*pi*X/Lx) * np.sin(3*pi*Y/Ly) / np.cosh(X-Lx/2)**2 / np.cosh(Y-Ly/2)**2
p0 = 0*q0
Q, P, T = wave2D(q0, p0, x, y, mask = c_mask, c = 1, t0 = 0, tf = 20, g = .4, boundary = 'fixed')
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
anim = animate2d(Q, x, y, T, fps = 24, frames = 144)
display(anim)