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
semh = SEMhat(order)
N = 100
n_dofs = (order+1)*N
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])
EtoV[-1,-1] = EtoV[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
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 [ ]:
# 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 [ ]:
# Time stepping
def g(u):
f = calc_flux(u)
c = np.max(calc_eig(u))
flux = AVG.dot(f)+c/2.0*JUMP.dot(u)
return Binv.dot(-S.dot(f)+SD.dot(FD.dot(f)-flux))
In [ ]:
# Integrate with SSP-RK3
u = u0.copy()
# 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)
v2 = .25*(3*u+v1+dt*g(v1))
u = (u+2*v2+2*dt*g(v2))/3.0
# 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
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)
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 [ ]: