In [ ]:
%matplotlib inline
%load_ext autoreload
%autoreload 2
In [ ]:
import numpy as np
from numpy import newaxis
import scipy.sparse as sps
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
In [ ]:
from fem import *
eval_phi1d = eval_lagrange_d0
In [ ]:
from poly import eval_P
from utils import minmod
from dg import *
In [ ]:
N = 1
periodic = False
#semh = SEMhat(order)
K = 250
n_dofs = (N+1)*K
do_limiter = True
do_filter = False
In [ ]:
mesh = Mesh(K, N)
In [ ]:
# Euler -- Sod's Shock Tube
use_filter = False
eqnset = EqnSetEuler()
gamma = eqnset.gamma
# Set up classical Sod's shock tube
u0 = np.zeros((n_dofs, 3))
u0[:,0] = 1.0
u0[:,0][mesh.dof_phys>=0.5] = 0.125
u0[:,1] = 0.0
u0[:,2] = 1.0
u0[:,2][mesh.dof_phys>=0.5] = 0.1
u0[:,2] /= (gamma-1)
ue = u0.copy()
np.max(eqnset.calc_eig(u0))
In [ ]:
# # Isothermal Euler -- Gaussian bump
# a = 1.0
# def calc_flux(u):
# f = np.zeros_like(u)
# f[:,0] = u[:,1]
# f[:,1] = u[:,1]**2/u[:,0]+a*a*u[:,0]
# f[:,2] = 0.0
# return f
# def calc_eig(u):
# return a+np.abs(u[:,1]/u[:,0])
# u0 = np.zeros((n_dofs, 3))
# aa = 100.*(dof_phys-L/2.0)**2
# u0[:,0] = 2.0+.2*np.exp(-aa)
# ue = u0.copy()
# np.max(calc_eig(u0))
In [ ]:
# # Shallow water -- Gaussian bump
# a = 1.0
# def calc_flux(u):
# f = np.zeros_like(u)
# f[:,0] = u[:,1]
# f[:,1] = u[:,1]**2/u[:,0]+.5*a*u[:,0]**2
# f[:,2] = 0.0
# return f
# def calc_eig(u):
# return u[:,1]/u[:,0]+np.sqrt(a*u[:,0])
# u0 = np.zeros((n_dofs, 3))
# aa = 100.*(dof_phys-L/2.0)**2
# u0[:,0] = 1.0+.5*np.exp(-aa)
# ue = u0.copy()
# np.max(calc_eig(u0))
In [ ]:
# # Burger -- Simple Shock
# def calc_flux(u):
# return u*u/2.0
# def calc_eig(u):
# return np.abs(u[:,0])
# ul = 2.0
# ur = 1.0
# s = (calc_flux(ur)-calc_flux(ul))/(ur-ul)
# u0 = np.ones((n_dofs, 3))*ur
# u0[dof_phys<=L/2.0,:] = ul
# # ue = np.ones_like(u0)*ur
# # ue[dof_phys<=L/2.0+s*dt*nt,:] = ul
# ue = u0.copy()
In [ ]:
# # Linear Advection -- Gaussian bump
# a = -2.0
# def calc_flux(u):
# return a*u
# def calc_eig(u):
# return np.abs(a)
# u0 = np.zeros((n_dofs, 3))
# aa = 100.*(dof_phys-L/2.0)**2
# u0[:,0] = np.exp(-aa)
# ue = u0.copy()
# np.max(calc_eig(u0))
In [ ]:
# Slope limiter
if do_limiter:
limiter = LimiterMUSCL(mesh)
else:
limiter = LimiterNULL(mesh)
def set_bc(u):
if not periodic:
u[[0,-1],:] = u0[[0,-1],:]
In [ ]:
# Filter
ss = 10.0
aa = 36.0
eta = np.arange(N+1.0)/N
SIG = sps.dia_matrix((np.exp(-aa*eta**ss),0),
shape=mesh.V.shape)
FILTER = sps.kron(sps.eye(N), mesh.V.dot(SIG.dot(mesh.Vinv)))
In [ ]:
Tmax = 0.2
CFL = 0.1
In [ ]:
# Time stepping
f0 = eqnset.calc_flux(u0[[0,-1],:])
def g(u):
f = eqnset.calc_flux(u)
c = np.max(eqnset.calc_eig(u))
AVG, JUMP = mesh.AVG, mesh.JUMP
Binv, C = mesh.Binv, mesh.C
SD, FD = mesh.SD, mesh.FD
flux = AVG.dot(f)+c/2.0*JUMP.dot(u)
if do_filter:
flux = FILTER.dot(flux)
if not periodic:
flux[[0,-1],:] = f0
return Binv.dot(-C.dot(f)+SD.dot(FD.dot(f)-flux))
In [ ]:
# Integrate with SSP-RK3
dx = mesh.dx
u = u0.copy()
f0 = eqnset.calc_flux(u0[[0,-1],:])
lambda_max = np.max(eqnset.calc_eig(u0))
meig = []
t = 0.0
while t<Tmax:
dt = CFL*dx/lambda_max
if t+dt>Tmax:
dt = Tmax-t
t += dt
v1 = u+dt*g(u)
set_bc(v1)
limiter.apply_limiter(v1)
v2 = .25*(3*u+v1+dt*g(v1))
set_bc(v2)
limiter.apply_limiter(v2)
u = (u+2*v2+2*dt*g(v2))/3.0
set_bc(u)
limiter.apply_limiter(u)
meig += [np.max(eqnset.calc_eig(u))]
meig = np.array(meig)
In [ ]:
k = -1 if periodic else n_dofs
dof_phys = mesh.dof_phys
plt.figure()
plt.plot(dof_phys[:k], u[:,0][:k])
plt.plot(dof_phys[:k], ue[:,0][:k], 'g--')
plt.ylabel('$\\rho$', size=16)
plt.figure()
plt.plot(dof_phys[:k], (u[:,1]/u[:,0])[:k])
plt.plot(dof_phys[:k], ue[:,1][:k], 'g--')
plt.ylabel('$u$', size=16)
# plt.figure()
# plt.plot(dof_phys[:k], u[:,2][:k])
# plt.plot(dof_phys[:k], ue[:,2][:k], 'g--')
# plt.ylabel('$E$', size=16)
In [ ]:
np.sum(CFL*dx/meig)*np.max(meig)
In [ ]:
plt.plot(CFL*dx/meig)
In [ ]: