In [1]:
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy.fftpack import fft, ifft,fft2,fftshift
import datetime
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [2]:
print('Last modified =', "{0:%Y-%b-%d %H:%M:%S.%f}".format(datetime.datetime.now()))
The Linear theory¶
Considering a multi-component plasma with different temperatures between $z$-direction($T_{||}$) and $x, y$-direction($T_\perp$). Writing the temperature ratio for $s$ species as $\alpha_s^2=T_{s,||}/T_{s,\perp}$ and the $z$-direction thermal velocity $v_{th,s}$, a distribution function is written as follows, $${f_s}\left( {{v_s},{v_y},{v_z}} \right) = \frac{{{n_s}}}{{{\pi ^{3/2}}{\alpha _s}v_{th,s}^3}}\exp \left[ { - \frac{{v_x^2 + v_y^2}}{{v_{th,s}^2}} - \frac{{v_z^2}}{{\alpha _s^2v_{th,s}^2}}} \right]$$. Here, $n_s$ is the number density.
When there is no background magnetic field, the dispersion relation of an electromagnetic mode with wavenumber vectors perpendicular to $z$-axis is given by,
$${\omega ^2} - {\left( {kc} \right)^2} + \sum\limits_s {\omega _{ps}^2\left[ {{A_s} + \left( {{A_s} + 1} \right){\zeta _s}Z\left( {{\zeta _s}} \right)} \right]} = 0.$$
Here, $\omega_{ps}=\sqrt{\frac{4\pi n_sq_s^2}{m_s}}$, $\zeta_s=\zeta_s(\omega,k)=\frac{\omega}{kv_{th,s}}$ and $Z(\zeta)$ is the plasma dispertion function,
$$Z(\zeta ) \equiv \frac{1}{\pi }\int_{ - \infty }^\infty {\frac{{{{\text{e}}^{ - {z^2}}}}}{{z - \zeta }}dz}$$.
Moreover, $A_s$ is a parameter indicating the magnitude of temperature anisotropy,
$${A_s} = \frac{{{T_{s,||}}}}{{{T_{s, \bot }}}} - 1 = \alpha _s^2 - 1$$.
Simulation results¶
Parameters¶
- electron-positron plasma
- $v_{th,e^-}=v_{th,e^+}=0.1c$
- $\alpha=5$
- $N_p=50$
- $(L_x, L_y)=(40\lambda_{De}, 40\lambda_{De})$
- periodic boundary
- $\Delta x=c=m_e=1$
- $\Delta t=0.4\Delta x$
In [3]:
#parameters
nx=256
ny=16
nrank=16
nt=2000
lx=40
ly=40
intvl=5
nd=nt//intvl #the number of data
rmass=1
dt=0.4*0.156
dir='data' #the directory location of the data
Snapshot of the electromagnetic fields¶
In [4]:
%%time
it=500
nmv=['bx','by','bz','ex','ey','ez']
bx,by,bz,ex,ey,ez=np.zeros((6,ny*nrank,nx))
EBs=[bx,by,bz,ex,ey,ez]
for i,EB in enumerate(EBs):
for irank in range(nrank):
f=open('%s/%s%06d_rank=%03d.dat' %(dir,nmv[i],it,irank),'rb')
d=np.fromfile(f,dtype='float64').reshape(-1,nx+2)
EB[ny*irank:ny*(irank+1),:]=d[0:ny,0:nx]
In [5]:
%%time
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(16,8))
plt.subplots_adjust(hspace=0.5,wspace=0.2)
im1=ax1.imshow(bx,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im1,ax=ax1);ax1.set_title(r'$B_x$')
im2=ax2.imshow(by,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im2,ax=ax2);ax2.set_title(r'$B_y$')
im3=ax3.imshow(bz,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im3,ax=ax3);ax3.set_title(r'$B_z$')
im4=ax4.imshow(ex,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im4,ax=ax4);ax4.set_title(r'$E_x$')
im5=ax5.imshow(ey,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im5,ax=ax5);ax5.set_title(r'$E_y$')
im6=ax6.imshow(ez,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im6,ax=ax6);ax6.set_title(r'$E_z$')
ax1.set_xlabel(r'$x/\lambda_D$');ax2.set_xlabel(r'$x/\lambda_D$');ax3.set_xlabel(r'$x/\lambda_D$');ax4.set_xlabel(r'$x/\lambda_D$');ax5.set_xlabel(r'$x/\lambda_D$');ax6.set_xlabel(r'$x/\lambda_D$')
ax1.set_ylabel(r'$y/\lambda_D$');ax2.set_ylabel(r'$y/\lambda_D$');ax3.set_ylabel(r'$y/\lambda_D$');ax4.set_ylabel(r'$y/\lambda_D$');ax5.set_ylabel(r'$y/\lambda_D$');ax6.set_ylabel(r'$y/\lambda_D$')
plt.show()
In [6]:
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(15,5))
plt.subplots_adjust(hspace=0.4)
im1=ax1.plot(np.mean(bx,axis=0));ax1.set_title('Bx')
im2=ax2.plot(np.mean(by,axis=0));ax2.set_title('By')
im3=ax3.plot(np.mean(bz,axis=0));ax3.set_title('Bz')
im4=ax4.plot(np.mean(ex,axis=0));ax4.set_title('Ex')
im5=ax5.plot(np.mean(ey,axis=0));ax5.set_title('Ey')
im6=ax6.plot(np.mean(ez,axis=0));ax6.set_title('Ez')
ax4.set_ylim(-1,1)
plt.show()
In [7]:
plt.imshow(np.sqrt(bx**2+by**2),origin='lower',aspect='auto',cmap='gist_heat')
plt.colorbar(label=r'$|B|$');
plt.title(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
plt.show()
Time evolution of the magnetic fields¶
In [8]:
%%time
tbz=np.zeros(nd)
for it in range(nd):
iit=intvl*(it+1)
for irank in range(nrank):
f1=open('%s/bx%06d_rank=%03d.dat' %(dir,iit,irank),'rb')
f2=open('%s/by%06d_rank=%03d.dat' %(dir,iit,irank),'rb')
d1=np.fromfile(f1,dtype='float64').reshape(-1,nx+2)
d2=np.fromfile(f2,dtype='float64').reshape(-1,nx+2)
tbz[it]=np.mean(np.sqrt(d1**2+d2**2))
plt.plot(tbz,'k')
plt.show()
Animation¶
In [ ]:
!mkdir figures
In [ ]:
%%time
nmv=['bx','by','bz']
B=np.zeros((3,ny*nrank,nx))
for iit in range(20):
it=(iit+1)*intvl
for i in range(3):
for irank in range(nrank):
f=open('%s/%s%06d_rank=%03d.dat' %(dir,nmv[i],it,irank),'rb')
d=np.fromfile(f,dtype='float64').reshape(-1,nx+2)
B[i,ny*irank:ny*(irank+1),:]=d[0:ny,0:nx]
fig= plt.subplots(figsize=(8,8))
#plt.subplots_adjust(hspace=0.5,wspace=0.3)
#fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
plt.imshow(np.sqrt(B[0,:,:]**2+B[1,:,:]**2+B[2,:,:]**2),aspect='auto',origin='lower',cmap='gist_heat',extent=[0,lx,0,ly])
plt.colorbar(label=r'$|B|$');
plt.title(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
plt.xlabel(r'$x/\lambda_D$')
plt.ylabel(r'$y/\lambda_D$')
plt.savefig('figures/B%03d.png' %(iit), bbox_inches='tight')
#print(iit)
plt.close()
In [ ]:
!convert -delay 5 -loop 0 figures/B*.png Bmag.gif
Animation is here!¶
fio.f90¶
Part of fio.f90 is modified as follows
write(*,*) 'nxs=',nxs,'nxe=',nxe,'nys=',nys,'nye=',nye
write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'bx',it2,'_rank=',nrank,'.dat'
open(110+nrank,file=fname,form='unformatted',access="stream")
do i=nys-1,nye+1
write(110+nrank) (uf(1,j,i), j =nxs-1,nxe+1)
end do
close(110+nrank)
write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'by',it2,'_rank=',nrank,'.dat'
open(110+nrank,file=fname,form='unformatted',access="stream")
do i=nys-1,nye+1
write(110+nrank) (uf(2,j,i), j=nxs-1,nxe+1)
end do
close(110+nrank)
write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'bz',it2,'_rank=',nrank,'.dat'
open(110+nrank,file=fname,form='unformatted',access="stream")
do i=nys-1,nye+1
write(110+nrank) (uf(3,j,i), j=nxs-1,nxe+1)
end do
close(110+nrank)
write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'ex',it2,'_rank=',nrank,'.dat'
open(110+nrank,file=fname,form='unformatted',access="stream")
do i=nys-1,nye+1
write(110+nrank) (uf(4,j,i), j=nxs-1,nxe+1)
end do
close(110+nrank)
write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'ey',it2,'_rank=',nrank,'.dat'
open(110+nrank,file=fname,form='unformatted',access="stream")
do i=nys-1,nye+1
write(110+nrank) (uf(5,j,i), j=nxs-1,nxe+1)
end do
close(110+nrank)
write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'ez',it2,'_rank=',nrank,'.dat'
open(110+nrank,file=fname,form='unformatted',access="stream")
do i=nys-1,nye+1
write(110+nrank) (uf(6,j,i), j=nxs-1,nxe+1)
end do
close(110+nrank)
if(mod(it2,5000)==0) then
write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'up1',it2,'_rank=',nrank,'.dat'
open(110+nrank,file=fname,form='unformatted',access="stream")
do j=nys,nye
do ii=1,np2(j,1)
write(110+nrank) (up(i,ii,j,1), i=1,5)
end do
end do
close(110+nrank)
write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'up2',it2,'_rank=',nrank,'.dat'
open(110+nrank,file=fname,form='unformatted',access="stream")
do j=nys,nye
do ii=1,np2(j,2)
write(110+nrank) (up(i,ii,j,2), i=1,5)
end do
end do
close(110+nrank)
endif