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.

In [1]:
%%bash
cat input.f90

module input
!### definition by a user ###!
integer, parameter :: ppc=500
integer, parameter :: nx=512
integer, parameter :: nsp=2
integer, parameter :: nt=8192
integer, parameter :: isav=8
integer, parameter :: nc=10
real(8), parameter :: pi=3.141592653589793d0
real(8), parameter :: Te=0.01d0
real(8), parameter :: vt=0.1d0
real(8), parameter :: dx=0.5d0
real(8), parameter :: dt=0.01d0
real(8), parameter :: b0=1.d0
real(8), parameter :: th=0
real(8), parameter :: eta=1.d-4
real(8), parameter :: q(nsp)    =(/1.d0, 1.d0/)
real(8), parameter :: rm(nsp)   =(/1.d0, 1.d0/)
real(8), parameter :: qomdt(nsp)=(/0.5*dt*q(1)/rm(1), 0.5*dt*q(2)/rm(2)/)
real(8), parameter :: dn(nsp)   =(/0.80d0/ppc, 0.2d0/ppc/)
real(8), parameter :: vx0(nsp)  =(/ -0.75d0, 3.d0/)

!### don't change if you have no reason! ###!
integer, parameter :: npmax=ppc*nx
integer, dimension(nsp) :: np
integer, dimension(nx) :: jj, jp, jm
real(8), parameter :: lx=dx*nx



end module input

Resonant Ion Beam Instability

Parameters

  • $n_x=512$
  • $ppc=500$
  • $n_t=8192$
  • $d_x=0.5$
  • $d_t=0.01$
  • $n_{CL}=10$ (sub-cycle time steps)
  • $\beta_e=0.01$
  • $\beta_i=1$
  • $n_b=0.2$
  • $v_b=3$
  • $\theta_{B}=0^\circ$
In [36]:
%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,ifft2,fftshift,ifftshift
import time
from IPython import display
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [32]:
nx=512
nt=1024
isav=8
dx=0.5
dt=0.025*isav
lx=dx*nx
lt=dt*nt
In [25]:
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)

Snapshots

In [26]:
plt.figure(figsize=(12,12))
plt.subplot(221);plt.imshow(nihst,origin='lower',cmap='jet',aspect='auto',extent=(0,nx*dx,0,nt*dt))
plt.title(r'$\rho_i$');plt.xlabel(r'$x/(c/\omega_{pi})$');plt.ylabel(r'$t\Omega_{ci}$')
plt.subplot(222);plt.imshow(vixhst,origin='lower',cmap='jet',aspect='auto',extent=(0,nx*dx,0,nt*dt))
plt.title(r'$V_i$');plt.xlabel(r'$x/(c/\omega_{pi})$');plt.ylabel(r'$t\Omega_{ci}$');#plt.colorbar()
plt.subplot(223);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(224);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.show()

Time Evolution

In [49]:
tt=dt*np.arange(nt)
bhst=np.sqrt(byhst**2+bzhst**2)
plt.plot(tt,np.mean(exhst**2,axis=1),label=(r'$E_x^2$'))
plt.plot(tt,np.mean(eyhst**2,axis=1),label=(r'$E_y^2$'))
plt.plot(tt,np.mean(ezhst**2,axis=1),label=(r'$E_z^2$'))
plt.plot(tt,np.mean(byhst**2,axis=1),label=(r'$B_y^2$'))
plt.plot(tt,np.mean(bzhst**2,axis=1),label=(r'$B_z^2$'))
plt.legend()
plt.xlabel(r'$t\Omega_{ci}$');plt.ylabel(r'$|B|$');
plt.yscale('log')
plt.show()
In [41]:
%%time
plt.figure(figsize=(12,8))
fftbz=fftshift(fft2(bzhst[:,:]))
bzpls=np.zeros((nt,nx), dtype=np.complex)
bzmin=np.zeros((nt,nx), dtype=np.complex)
plt.subplot(131)
plt.imshow(bzhst,aspect='auto',origin='lower',extent=[0,lx,0,lt],clim=(-2,2),cmap='bwr')
plt.colorbar()
plt.title(r'$B_z^{raw}$')
plt.xlabel(r'$x/(c/\omega_{e})$')
plt.ylabel(r'$t\Omega_{ci}$')

plt.subplot(132)
bzpls[0:nt//2,nx//2:nx]=fftbz[0:nt//2,nx//2:nx]
bzpls[nt//2:nt,0:nx//2]=fftbz[nt//2:nt,0:nx//2]
bzpls=ifft2(ifftshift(bzpls))
plt.imshow(np.real(bzpls),aspect='auto',origin='lower',extent=[0,lx,0,lt],clim=(-1.5,1.5),cmap='bwr')
plt.colorbar()
plt.title(r'$B_z^+$')

plt.subplot(133)
bzmin[0:nt//2,0:nx//2]  =fftbz[0:nt//2,0:nx//2]
bzmin[nt//2:nt,nx//2:nx]=fftbz[nt//2:nt,nx//2:nx]
bzmin=ifft2(ifftshift(bzmin))
plt.imshow(np.real(bzmin),aspect='auto',origin='lower',extent=[0,lx,0,lt],clim=(-1,1),cmap='bwr')
plt.colorbar()
plt.title(r'$B_z^-$')

plt.show()

Wall time: 2.65 s