In [ ]:
%matplotlib inline
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
In [ ]:
order = 1
periodic = True
semh = SEMhat(order)
N = 100
n_dofs = (order+1)*N
do_filter = False
In [ ]:
L = 1.0
vertices = np.linspace(0, L, N+1)
EtoV = np.zeros((N, 2), dtype=np.int)
EtoV[:,0] = np.arange(N)
EtoV[:,1] = np.arange(N)+1
topo = Interval()
xq = topo.ref_to_phys(vertices[EtoV], semh.xgll)
jacb_det = topo.calc_jacb(vertices[EtoV])[0]
dx = np.min(xq[0,1:]-xq[0,:-1])
if periodic:
EtoV[-1,-1] = EtoV[0,0]
# Restriction operator
# Removes contribution to dirch points from flux operators
R = sps.eye(n_dofs).tocsr()
if not periodic:
R[0,0] = R[-1,-1] = 0.0
In [ ]:
# Make elem to dof map
EtoD = np.arange(N*(order+1))
EtoD = EtoD.reshape((N, -1))
dof_phys = xq.ravel()
# Averaging operator
rows = EtoD[:,[0,-1]].ravel()
cols = EtoV.ravel()
vals = np.ones_like(cols)
FtoD = sps.coo_matrix((vals, (rows, cols))).tocsr()
AVG = FtoD.dot(FtoD.T)/2.0
# Extract face DOFS
vals = np.ones(len(rows))
FD = sps.coo_matrix((vals, (rows, rows))).tocsr()
# Set face signs
vals[::2] = -1
SD = sps.coo_matrix((vals, (rows, rows))).tocsr()
# Jump operator
JUMP = FtoD.dot(SD.dot(FtoD).T)
In [ ]:
# Build Advection operator
S = sps.kron(sps.eye(N), semh.Ch).tocsr()
# Differentiation matrix
Dr = sps.kron(sps.eye(N), semh.Dh)/jacb_det
Dr = Dr.tocsr()
# Build full elemental mass matrix
x, w = topo.get_quadrature(order+1)
P = eval_phi1d(semh.xgll, x).T
G = sps.dia_matrix((w, 0), shape=(len(x), len(x)))
Bf = P.T.dot(G.dot(P))*jacb_det
Bfinv = np.linalg.inv(Bf)
B = sps.kron(sps.eye(N), Bf)
# Using trick from book
V = eval_P(order, semh.xgll).T
Vinv = np.linalg.inv(V)
Minv = V.dot(V.T)/jacb_det
Binv = sps.kron(sps.eye(N), Minv).tocsr()
In [ ]:
# Build Poisson matrix
tau = 1.0
FLUXU = AVG
Q = Dr-Binv.dot(SD.dot(FD-FLUXU))
FLUXQ = AVG.dot(Q)-tau*JUMP
A = S.dot(Q)-SD.dot(FD.dot(Q)-FLUXQ)
A = -A
In [ ]:
# # Euler -- Sod's Shock Tube
# use_filter = False
# gamma = 7.0/5.0
# def calc_flux(u):
# f = np.zeros_like(u)
# p = (gamma-1)*(u[:,2]-.5*u[:,1]**2/u[:,0])
# f[:,0] = u[:,1]
# f[:,1] = u[:,1]**2/u[:,0]+p
# f[:,2] = (u[:,2]+p)*u[:,1]/u[:,0]
# return f
# def calc_p(u):
# p = (gamma-1)*(u[:,2]-.5*u[:,1]**2/u[:,0])
# return p
# def calc_eig(u):
# p = calc_p(u)
# eig = np.abs(u[:,1]/u[:,0])
# eig += np.sqrt(gamma*p/u[:,0])
# return eig
# # Set up classical Sod's shock tube
# u0 = np.zeros((n_dofs, 3))
# u0[:,0] = 1.0
# u0[:,0][dof_phys>=0.5] = 0.125
# u0[:,1] = 0.0
# u0[:,2] = 1.0
# u0[:,2][dof_phys>=0.5] = 0.1
# u0[:,2] /= (gamma-1)
# ue = u0.copy()
# np.max(calc_eig(u0))
In [ ]:
# Isothermal Euler -- Gaussian bump
a = 0.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])
return 1.0
u0 = np.zeros((n_dofs, 3))
alpha = 0.01
u0[:,0] = 2.0+alpha*np.sin(2*np.pi*dof_phys)
wp = np.sqrt(np.mean(u0[:,0]))
ue = u0.copy()
np.max(calc_eig(u0))
In [ ]:
def phie(t):
return alpha*np.sin(2*np.pi*dof_phys)*np.cos(wp*t)/((2*np.pi)**2)
def Ee(t):
return -alpha*np.cos(2*np.pi*dof_phys)*np.cos(wp*t)/(2*np.pi)
In [ ]:
def set_bc(u):
if not periodic:
u[[0,-1],:] = u0[[0,-1],:]
def apply_limiter(u):
pass
In [ ]:
# Filter
ss = 16.0
aa = 36.0
eta = np.arange(order+1.0)/order
SIG = sps.dia_matrix((np.exp(-aa*eta**ss),0),
shape=V.shape)
FILTER = sps.kron(sps.eye(N), V.dot(SIG.dot(Vinv)))
In [ ]:
# Time stepping
f0 = calc_flux(u0[[0,-1],:])
def g(u):
f = calc_flux(u)
c = np.max(calc_eig(u))
flux = AVG.dot(f)+c/2.0*JUMP.dot(u)
if not periodic:
flux[[0,-1],:] = f0
return Binv.dot(-S.dot(f)+SD.dot(FD.dot(f)-flux))
In [ ]:
# Integrate with SSP-RK3
u = u0.copy()
f0 = calc_flux(u0[[0,-1],:])
lambda_max = np.max(calc_eig(u0))
# One complete oscillation
Tfinal = 2*np.pi/wp
dt = 0.001
nt = int(Tfinal/dt)
dt = Tfinal/nt
for k in range(nt):
# Step Euler equations
v1 = u+dt*g(u)
set_bc(v1)
apply_limiter(v1)
v2 = .25*(3*u+v1+dt*g(v1))
set_bc(v2)
apply_limiter(v2)
u = (u+2*v2+2*dt*g(v2))/3.0
set_bc(u)
apply_limiter(u)
# Step field equations
phi = sps.linalg.spsolve(A, B.dot(u[:,0]-np.mean(u[:,0])))
E = -Dr.dot(phi)
u[:,1] += dt*u[:,0]*E
phif = sps.linalg.spsolve(A, B.dot(u[:,0]-np.mean(u[:,0])))
Ef = -Dr.dot(phif)
nt, dt*nt
In [ ]:
k = -1 if periodic else n_dofs
plt.figure()
plt.plot(dof_phys[:k], ue[:,0][:k], 'g--',
label="$t=0$",
linewidth=2)
plt.plot(dof_phys[:k], u[:,0][:k],
label="$t=2\pi/\omega_p$")
plt.ylabel('$\\rho$', size=16)
plt.xlabel("$x$", size=16)
#plt.legend(loc='upper right', fontsize=16)
plt.title("One complete Langmuir Oscillation", size=16)
plt.savefig("langmuir.pdf")
In [ ]:
ps = lambda x:x.reshape((-1,order+1)).T
plt.figure()
plt.plot(ps(dof_phys), ps(u[:,1]/u[:,0]), 'b')
plt.plot(ps(dof_phys), ps(ue[:,1]), '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 [ ]:
plt.plot(ps(dof_phys), ps(phif-np.mean(phif)), 'b');
plt.plot(dof_phys, phie(nt*dt), 'g--',
linewidth=2)
In [ ]:
plt.plot(ps(dof_phys), ps(Ef), 'b');
plt.plot(dof_phys, Ee(nt*dt), 'g--',
linewidth=2)
In [ ]: