

In [1]:

    
%matplotlib inline
%load_ext autoreload
%autoreload 2



In [2]:

    
import numpy as np
from pic1d2v import *
from picwave import *
import matplotlib.pyplot as plt
from matplotlib import animation
from collections import namedtuple
norm = lambda x: np.max(np.abs(x))

Landau Oscillations

Validate pic1d2v by measuring numerical angular frequency of oscillations resulting from perturbed particle density, and comparing with the expected numerical dispersion relation



In [3]:

    
mode_all = 1
pw = PicWave(mode=mode_all)
pw.run()

print 'wp:', pw.wp, '\tq: ', pw.q, '\tm: ', pw.m









    



wp: 1.0 	q:  0.0245436926062 	m:  0.0245436926062

Error in each particles SHM (w.r.t dispersion free exact solution)



In [4]:

    
xt     = pw.delta_x()
t_vals = np.linspace(0,pw.dt*pw.nt, pw.nt+1)
exact  = lambda i: np.cos(pw.wp*t_vals)*xt[0, i]
errors = [norm(xt[:,i]-exact(i)) 
          for i in range(pw.N_e)]

plt.ylabel("Error"); plt.xlabel("Particle")
plt.plot(errors);



In [5]:

    
# Animation of particle motion over time
if False:
    fig = plt.figure()
    p1 = fig.add_subplot(1, 1, 1)
    p1.set_xlabel('x'); p1.set_ylabel('$v_x$')
    imgs = []   
    s = 1
    x_vals = np.linspace(0, pw.L, pw.nx+1)[:-1]
    for i in range(pw.nt+1):
        imgs.append((p1.scatter(pw.xp[i,:pw.N_e], pw.vx[i,:pw.N_e], 
                                color='b', s=s),
                     p1.plot(x_vals, pw.phi[i], color='g')[0],
                    ))


    im_ani = animation.ArtistAnimation(fig, imgs, 
                                       interval=50, repeat_delay=3000,
                                       blit=True)
    im_ani.save("1d2v-wave.mp4", writer="mencoder")



In [6]:

    
# Animation of k vals over time
pw2 = PicWave(mode=mode_all)
pw2.run()
k_vals = np.fft.rfft(pw2.delta_x(), axis=1)

if False:
    fig = plt.figure()
    p1 = fig.add_subplot(1, 1, 1)
    p1.set_xlabel('x'); p1.set_ylabel('$v_x$')
    imgs = []   
    for i in range(pw2.nt+1):
        imgs.append((p1.plot(np.abs(k_vals[i]), color='b')[0],))


    im_ani = animation.ArtistAnimation(fig, imgs, 
                                       interval=50, repeat_delay=3000,
                                       blit=True)
    im_ani.save("1d2v-wavek.mp4", writer="mencoder")

$\omega(k)$ measured vs. expected as a function of $k$

Varying $k$



In [7]:

    
modes = range(1, 50)
t_vals = np.linspace(0, pw.dt*pw.nt, pw.nt+1)
points = np.arange(1, 6)
k_vals = []
wk_measured = []
wk_expected = []
for mode in modes:
    pwm = PicWave(mode=mode)
    pwm.run()
    rhok = np.fft.rfft(pwm.rho, axis=1)
    phik = np.fft.rfft(pwm.phi, axis=1)
    ESEk = np.abs(rhok[:, mode]*phik[:, mode].conj())
    ESEk2 = ESEk/ESEk[0] # Scale by BC (A(0))
    
    omegak = np.mean(np.arccos(2*ESEk2[points]-1)/t_vals[points]/2.)
    wk_measured.append(omegak)
    wk_expected.append(pwm.wp*np.cos(pwm.k*pwm.dx/2))
    k_vals.append(pwm.k)

plt.plot(k_vals, wk_measured, label="$\omega_k$ measured")
plt.plot(k_vals, wk_expected, label="$\omega_k$ expected")
plt.xlabel("$k$"); plt.ylabel("$\omega_k$")
plt.legend();

Metrics for determining model parameters

I considered the following model parameters

$kdx$
particles per cell
Time step $dt$

$\omega(k)$ measured vs. expected as a function of $k dx$

Holding $k$ constant and varying $dx$ (but not $L$)



In [8]:

    
nx_vals  = [2**k for k in range(6, 12+1)]
mode     = 20
# Rough number of particles per cell
ppc_vals = [.5, 1, 2, 3]
t_vals = np.linspace(0, pw.dt*pw.nt, pw.nt+1)
points = np.arange(1, 6)
for ppc in ppc_vals:
    kdx_vals    = []
    wk_measured = []
    wk_expected = []
    for nx in nx_vals:
        pwm = PicWave(mode=mode, nx=nx, N_e=int(nx*ppc))
        pwm.run()
        rhok = np.fft.rfft(pwm.rho, axis=1)
        phik = np.fft.rfft(pwm.phi, axis=1)
        ESEk = np.abs(rhok[:, mode]*phik[:, mode].conj())
        ESEk2 = ESEk/ESEk[0] # Scale by BC (A(0))

        omegak = np.mean(np.arccos(2*ESEk2[points]-1)/t_vals[points]/2.)
        wk_measured.append(omegak)
        wk_expected.append(pwm.wp*np.cos(pwm.k*pwm.dx/2))
        kdx_vals.append(pwm.k*pwm.dx)

    plt.plot(kdx_vals, wk_measured, label="ppc=%s" % ppc)

plt.plot(kdx_vals, wk_expected, 'bd', label="$\omega_k$ expected")
plt.xlabel("$k dx$"); plt.ylabel("$\omega_k$")
plt.legend(loc="lower left");

$\omega(k)$ as a function of k and particles per cell

ppc = average number of particles per cell with no density pertubation



In [9]:

    
modes = range(1, 30)
nx = 128
ppc_vals = [.5, 1, 2, 3]
t_vals = np.linspace(0, pw.dt*pw.nt, pw.nt+1)
points = np.arange(1, 6)
for ppc in ppc_vals:
    k_vals = []
    wk_measured = []
    wk_expected = []
    for mode in modes:
        pwm = PicWave(mode=mode, N_e=int(nx*ppc))
        pwm.run()
        rhok = np.fft.rfft(pwm.rho, axis=1)
        phik = np.fft.rfft(pwm.phi, axis=1)
        ESEk = np.abs(rhok[:, mode]*phik[:, mode].conj())
        ESEk2 = ESEk/ESEk[0] # Scale by BC (A(0))

        omegak = np.mean(np.arccos(2*ESEk2[points]-1)/t_vals[points]/2.)
        wk_measured.append(omegak)
        wk_expected.append(pwm.wp*np.cos(pwm.k*pwm.dx/2))
        k_vals.append(pwm.k)

    plt.plot(k_vals, wk_measured, label="ppc=%s" % ppc)

plt.plot(k_vals, wk_expected, 'bd', label="$\omega_k$ expected")
plt.xlabel("$k$"); plt.ylabel("$\omega_k$")
plt.legend(loc="lower left");

$\omega(k)$ as a function of k and time step $dt$

Expect no change until leap-frog becomes unstable



In [10]:

    
modes = range(1, 30)
nx = 128
T = 30.
nt_vals = [50, 100, 150, 200, 250]
dt_vals = [T/nt for nt in nt_vals]
points = np.arange(1, 6)
for nt, dt in zip(nt_vals, dt_vals):
    k_vals = []
    wk_measured = []
    wk_expected = []
    for mode in modes:
        pwm = PicWave(mode=mode, nt=nt, dt=dt, N_e=int(nx*2))
        pwm.run()
        rhok = np.fft.rfft(pwm.rho, axis=1)
        phik = np.fft.rfft(pwm.phi, axis=1)
        ESEk = np.abs(rhok[:, mode]*phik[:, mode].conj())
        ESEk2 = ESEk/ESEk[0] # Scale by BC (A(0))
        
        t_vals = np.linspace(0, pwm.dt*pwm.nt, pwm.nt+1)
        omegak = np.mean(np.arccos(2*ESEk2[points]-1)/t_vals[points]/2.)
        wk_measured.append(omegak)
        wk_expected.append(pwm.wp*np.cos(pwm.k*pwm.dx/2))
        k_vals.append(pwm.k)

    plt.plot(k_vals, wk_measured, label="dt=%.2f" % dt)

plt.plot(k_vals, wk_expected, 'bd', label="$\omega_k$ expected")
plt.xlabel("$k$"); plt.ylabel("$\omega_k$")
plt.legend(loc="lower left");



In [10]: