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 *

Build Mesh and Operators


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)

Problem Setup


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))

Slope Limiter


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],:]

Filter


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)))

Compute solution


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 [ ]: