In [403]:
import numpy as np
import matplotlib.pyplot as plt
In [424]:
def u0(x, flag=1):
if flag == 0:
return np.exp(-10*((x-0.5)**2))
if flag == 1:
return (0.25*np.sin(2*np.pi*x)+0.5*np.sin(4*np.pi*x)+0.75)
if flag == 2:
ee = 0.005
return (1+0.5*(1-np.tanh((x-0.3)/4/ee))-0.5*(1-np.tanh((x-0.1)/4/ee)))+ \
(-1+0.5*(1-np.tanh((x-0.7)/4/ee))-0.5*(1-np.tanh((x-0.9)/4/ee)))
def get_flux(u, flag=1):
u_jp1 = np.roll(u,-1)
u_jm1 = np.roll(u,1)
#if flag == 0:
# return (((u + u_jp1)/2)**2 - ((u + u_jm1)/2)**2)
if flag == 0:
return (u_jp1**2 + 2*(u_jp1 -u_jm1)*u-u_jm1**2)/4
elif flag == 1:
u_rb = 0.5*(u+u_jp1)
u_r = (u**2)*(u_rb>=0) + (u_jp1**2)*(u_rb<0)
return u_r - np.roll(u_r,1)
elif flag == 2:
u_jp2 = np.roll(u,-2)
u_rb = 0.5*(u+u_jp1)
u_r = ((6*u/8 + 3*u_jp1/8 - u_jm1/8)*(u_rb>=0) + \
(6*u_jp1/8 + 3*u/8 - u_jp2/8)*(u_rb<0))**2
return u_r - np.roll(u_r,1)
elif flag == 3:
u_l = ((u + u_jm1)/2)**2 + np.maximum(abs(u),abs(u_jm1))*(u_jm1 - u)/2
return -u_l + np.roll(u_l,-1)
def step(u,Δt, Δx, flag):
return u - Δt*get_flux(u,flag)/(2*Δx)
In [464]:
def plot(axs, init_cond, scheme, t_save=[0,0.25,0.5,1]):
Nx = 250
Nt = 10000
T = np.linspace(0,1,Nt+1)[1:]
Δt = T[1] - T[0]
x = np.linspace(0,1,Nx+1)[:-1]
Δx = x[1] - x[0]
saveCount = 1
u = u0(x,init_cond)
axs[0].plot(x,u)
for i,t in enumerate(T):
u = step(u, Δt, Δx, scheme)
if t>=t_save[saveCount]:
axs[saveCount].plot(x,u)
if scheme == 0:
axs[0].set_title("$t=0$")
axs[saveCount].set_title("$t={}$".format(t))
saveCount += 1
In [473]:
scheme_names = ["Central","Upwind","QUICK","Roe"]
init_cond_expr = ["$e^{-10(x-0.5)^{2}}$",
"$0.25\sin{(2\pi x)} + 0.5\sin{(4\pi x)} + 0.75$",
"Look-alike solitons"]
for init_cond,init_cond_name in enumerate(init_cond_expr):
fig, axs = plt.subplots(4,4,sharex=True,sharey=True, figsize=(8,6))
fig.suptitle("Initial condition: "+init_cond_name,fontsize=14)
u_init = u0(x,init_cond)
axs[0,0].set_ylim([np.min(u_init)-0.1,np.max(u_init)+0.1])
for ax in axs.flatten():
ax.xaxis.set_major_locator(plt.MaxNLocator(5))
ax.yaxis.set_major_locator(plt.MaxNLocator(5))
ax.grid(True)
for tick in ax.get_xticklabels():
tick.set_rotation(20)
for scheme, scheme_name in enumerate(scheme_names):
plot(axs[scheme], init_cond, scheme)
axs[scheme][0].set_ylabel(scheme_name, fontsize=14)
fig.savefig("q2_{}.pdf".format(init_cond))
