1D Hybrid Simulation

Basi Equations

\begin{align} \frac{{d{{\mathbf{x}}_i}}}{{dt}} &= {{\mathbf{v}}_i} \hfill \\ \frac{{d{{\mathbf{v}}_i}}}{{dt}} &= \frac{{{q_i}}}{{{m_i}}}\left( {{\mathbf{E}} + {{\mathbf{v}}_i} \times {\mathbf{B}}} \right) \hfill \\ \frac{{\partial {\mathbf{B}}}}{{\partial t}} &= - c\nabla \times {\mathbf{E}} \hfill \\ \nabla \times {\mathbf{B}} &= \frac{{4\pi }}{c}{\mathbf{J}} = \frac{{4\pi }}{c}{q_i}{n_i}\left( {{{\mathbf{V}}_i} - {{\mathbf{V}}_e}} \right) \hfill \\ e{n_e} &= {q_i}{n_i} \hfill \\ {n_e}{m_e}\frac{{d{{\mathbf{V}}_e}}}{{dt}} &= 0 = - \frac{{{{\mathbf{V}}_i} \times {\mathbf{B}}}}{c} - \frac{{\nabla {P_e}}}{{{q_i}{n_i}}} - \frac{{{\mathbf{B}} \times \left( {\nabla \times {\mathbf{B}}} \right)}}{{4\pi {q_i}{n_i}}} \hfill \\ \end{align}

We solve those equations with CAM_CL method [Matthews1994JCP].
In this simulation, $v_A$ (Alfv\'en speed) is the unit of speed, $c/\omega_{pi}=v_A/\Omega_{ci}$ (ion inertial length) is the unit of length, $\Omega_{ci}^{-1}$ (cyclotron time) is the unit of time.

Linear Waves

Parameters

  • $n_x=512$
  • $ppc=400$
  • $n_t=8192$
  • $d_x=0.5$
  • $d_t=0.1$
  • $n_{CL}=10$ (sub-cycle time steps)
  • $\beta_e=0.02$
  • $\beta_i=10^{-4}$
  • $\theta_{B}=0^\circ$
In [13]:
%config InlineBackend.figure_format = 'retina'
from numba import jit,f8,u8
import numpy as np
import math as mt
import matplotlib.pyplot as plt
import math as mt
from scipy.fftpack import fft, ifft,fft2,fftshift
import time
from IPython import display
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [18]:
nx=512
nt=1024
isav=8
dx=0.5
dt=0.1*isav
In [10]:
f1 =open('data/exhst.d'  ,'r')
f2 =open('data/eyhst.d'  ,'r')
f3 =open('data/ezhst.d'  ,'r')
f4 =open('data/bxhst.d'  ,'r')
f5 =open('data/byhst.d'  ,'r')
f6 =open('data/bzhst.d'  ,'r')
f7 =open('data/nihst.d'  ,'r')
f8 =open('data/vixhst.d' ,'r')
f9 =open('data/viyhst.d' ,'r')
f10=open('data/vizhst.d' ,'r')

exhst =np.fromfile(f1, dtype='float64',sep=" ").reshape(-1,nx)
eyhst =np.fromfile(f2, dtype='float64',sep=" ").reshape(-1,nx)
ezhst =np.fromfile(f3, dtype='float64',sep=" ").reshape(-1,nx)
bxhst =np.fromfile(f4, dtype='float64',sep=" ").reshape(-1,nx)
byhst =np.fromfile(f5, dtype='float64',sep=" ").reshape(-1,nx)
bzhst =np.fromfile(f6, dtype='float64',sep=" ").reshape(-1,nx)
nihst =np.fromfile(f7, dtype='float64',sep=" ").reshape(-1,nx)
vixhst=np.fromfile(f8, dtype='float64',sep=" ").reshape(-1,nx)
viyhst=np.fromfile(f9, dtype='float64',sep=" ").reshape(-1,nx)
vizhst=np.fromfile(f10,dtype='float64',sep=" ").reshape(-1,nx)
In [26]:
kmin=2*mt.pi/(dx*nx)*(-nx/2); kmax=2*mt.pi/(dx*nx)*(nx/2)
wmin=2*mt.pi/(dt*nt)*(-nt/2); wmax=2*mt.pi/(dt*nt)*(nt/2)
kaxis=np.linspace(kmin,kmax,nx); waxis=np.linspace(wmin,wmax,nt)

fftby=fftshift(fft2(exhst[:,::-1]))
fftbz=fftshift(fft2(bzhst[:,::-1]))

plt.figure(figsize=(12,12))
plt.subplot(221);plt.imshow(exhst,origin='lower',cmap='bwr',aspect='auto',extent=(0,nx*dx,0,nt*dt))
plt.title(r'$E_x$');plt.xlabel(r'$x/(c/\omega_{pi})$');plt.ylabel(r'$t\Omega_{ci}$')
plt.subplot(222);plt.imshow(byhst,origin='lower',cmap='bwr',aspect='auto',extent=(0,nx*dx,0,nt*dt))
plt.title(r'$B_z$');plt.xlabel(r'$x/(c/\omega_{pi})$');plt.ylabel(r'$t\Omega_{ci}$')

plt.subplot(223);plt.imshow(abs(fftby),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.axis([0,kmax/2,0,wmax/2])
plt.xlabel(r'$k(c/\omega_{pi})$');plt.ylabel(r'$\omega/\Omega_{ci}$')
plt.clim(0,1)
plt.subplot(224);plt.imshow(abs(fftbz),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.xlabel(r'$k(c/\omega_{pi})$');plt.ylabel(r'$\omega/\Omega_{ci}$')
plt.axis([kmin,kmax,0,wmax])
plt.clim(0,0.5)
plt.show()