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

Build Mesh and Operators


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

Problem Setup


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

Filter


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

Compute solution


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