Hall-MHD¶
Reference: Nariyuki2007NPG
This post shows Hall-MHD simulations and also linear waves and the modulational instability as an application.
Governing Equations¶
\begin{equation} \begin{gathered} {\partial _t}\rho + \nabla \cdot \left( {\rho {\mathbf{v}}} \right) = 0; \hfill \\ \rho {\partial _t}{\mathbf{v}} + \rho {\mathbf{v}} \cdot \nabla {\mathbf{v}} = - \nabla \left( {p + {{\left| {\mathbf{b}} \right|}^2}/2} \right) + \left( {{\mathbf{b}} \cdot \nabla } \right){\mathbf{b}}; \hfill \\ {\partial _t}{\mathbf{b}} = \nabla \times \left( {{\mathbf{v}} \times {\mathbf{b}}} \right) - \nabla \times \left( {{\rho ^{ - 1}}\left( {\nabla \times {\mathbf{b}}} \right) \times {\mathbf{b}}} \right); \hfill \\ \nabla \cdot {\mathbf{b}} = 0; \hfill \\ p = \beta \rho . \hfill \\ \end{gathered} \end{equation}1D system¶
We consider 1D case and those equations can be re-written as follows:
\begin{equation} \begin{gathered} {\partial _t} + {\partial _x}\left( {\rho u} \right) = 0; \hfill \\ {\partial _t}u + u{\partial _x}u = - \frac{1}{\rho }{\partial _x}\left( {p + {{\left| {\mathbf{b}} \right|}^2}/2} \right); \hfill \\ {\partial _t}v + u{\partial _x}v = \frac{{{1}}}{\rho }{\partial _x}{b_y}; \hfill \\ {\partial _t}w + u{\partial _x}w = \frac{{{1}}}{\rho }{\partial _x}{b_z}; \hfill \\ {\partial _t}{b_y} = - {\partial _x}\left( {u{b_y} - v - \frac{{{b_x}}}{\rho }{\partial _x}{b_z}} \right); \hfill \\ {\partial _t}{b_z} = - {\partial _x}\left( {u{b_z} - w + \frac{{{b_x}}}{\rho }{\partial _x}{b_y}} \right); \hfill \\ p = \beta \rho . \hfill \\ \end{gathered} \end{equation}We solve these equations numerically by means of the spectral method. They are normalized in a way that the units of time and space reciprocal of ion cyclotron frequency, and the ion intertial length, respectively.
In [1]:
%config InlineBackend.figure_format='retina'
from IPython.core.display import display, HTML, clear_output
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as sf
import matplotlib.animation as animation
from IPython.display import display,Markdown
plt.rcParams['font.size'] = 12
plt.rcParams['axes.linewidth'] = 1.5
In [2]:
def HallMHD(nx,nt,dx,dt,beta,bx,U,isav):
global jp1,jm1,kx,kd,kd2
kx =2*np.pi/(nx*dx)*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
kd =np.r_[np.ones(nx//3),np.zeros(nx//3+nx%3),np.ones(nx//3)] #for de-aliasing
hst =np.zeros((6,nt//isav,nx))
hst[:,0,:]=U[:,:]
Uf=sf.fft(U,axis=1)
mu =2e-2*dx/np.pi
eta=2e-2*dx/np.pi
for it in range(1,nt):
#---add viscocity and resistivity term---#
Uf[1:4,:]=np.exp(-kx**2*dt*mu )*Uf[1:4,:]
Uf[4:6,:]=np.exp(-kx**2*dt*eta)*Uf[4:6,:]
#---4th order Runge-Kutta scheme---#
LK1=flx(Uf )
LK2=flx(Uf+0.5*dt*LK1)
LK3=flx(Uf+0.5*dt*LK2)
LK4=flx(Uf+ dt*LK3)
Uf=Uf+dt*(LK1+2*LK2+2*LK3+LK4)/6
#---3rd order TVD Runge-Kutta scheme---#
#LK1=flx(Uf)
#Uf1=Uf+dt*LK1
#LK2=flx(Uf1)
#Uf2=3/4*Uf+1/4*(Uf1+dt*LK2)
#LK3=flx(Uf2)
#Uf =1/3*Uf+2/3*(Uf2+dt*LK3)
if(it%isav==0):
hst[:,it//isav,:]=np.real(sf.ifft(Uf,axis=1))
return locals()
def flx(Uf):
ff=np.zeros((6,nx),dtype='complex128')
rho=np.real(sf.ifft(Uf[0,:]*kd))
u =np.real(sf.ifft(Uf[1,:]*kd))
v =np.real(sf.ifft(Uf[2,:]*kd))
w =np.real(sf.ifft(Uf[3,:]*kd))
by =np.real(sf.ifft(Uf[4,:]*kd))
bz =np.real(sf.ifft(Uf[5,:]*kd))
ux =np.real(sf.ifft(1j*kx*Uf[1,:]*kd))
vx =np.real(sf.ifft(1j*kx*Uf[2,:]*kd))
wx =np.real(sf.ifft(1j*kx*Uf[3,:]*kd))
byx=np.real(sf.ifft(1j*kx*Uf[4,:]*kd))
bzx=np.real(sf.ifft(1j*kx*Uf[5,:]*kd))
ff[0,:]=-1j*kx*sf.fft(rho*u)
ff[1,:]=-sf.fft(u*ux+bx*by*byx/rho+bx*bz*bzx/rho)-1j*kx*sf.fft(np.log(rho))*beta
ff[2,:]=-sf.fft(u*vx-bx*byx/rho)
ff[3,:]=-sf.fft(u*wx-bx*bzx/rho)
ff[4,:]=-1j*kx*sf.fft(u*by-v*bx-bx*bzx/rho)
ff[5,:]=-1j*kx*sf.fft(u*bz-w*bx+bx*byx/rho)
return ff
Linear Waves¶
Linearizing and Fourier transforming, we get the dispersion realtion of them. \begin{equation} {k^2} = \frac{{{\omega ^2}}}{{1 \pm {\omega ^2}}} \end{equation}
In [3]:
%%time
nx=1024;nt=1024*40;dx=0.25;dt=0.1*dx;isav=10
U=np.zeros((6,nx))
x=np.arange(nx)*dx
U[0,:]=1
for ik in range(256):
ikk=ik-128
U[0,:]=U[0,:]+1e-4*np.random.randn()*np.sin(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())
U[4,:]=U[4,:]+1e-4*np.random.randn()*np.sin(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())#U[4,jm1]/4+U[4,:]/2+U[4,jp1]/4
U[5,:]=U[5,:]+1e-4*np.random.randn()*np.cos(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())#U[4,jm1]/4+U[4,:]/2+U[4,jp1]/4
bx=1; beta=0.5
data=HallMHD(nx,nt,dx,dt,beta,bx,U,isav)
locals().update(data)
Time-Space Evolution¶
In [4]:
plt.figure(figsize=(10,5))
plt.subplot(121);plt.imshow(hst[0,:,:],origin='lower',aspect='auto', extent=([0,nx*dx,0,nt*dt] ));plt.xlabel('x');plt.ylabel('t');plt.title(r'$\rho$')
plt.subplot(122);plt.imshow(hst[5,:,:],origin='lower',aspect='auto',cmap='bwr',extent=([0,nx*dx,0,nt*dt]),clim=(-4e-3,4e-3));plt.xlabel('x');plt.ylabel('t');plt.title(r'$b_z$')
plt.show()
Dispersion Relation¶
In [5]:
ddt=dt*isav
ddx=dx
nd=nt//isav
kmin=2*np.pi/(ddx*nx)*(-nx/2); kmax=2*np.pi/(ddx*nx)*(nx/2)
wmin=2*np.pi/(ddt*nd)*(-nd/2); wmax=2*np.pi/(ddt*nd)*(nd/2)
kaxis=np.linspace(kmin,kmax,nx); waxis=np.linspace(wmin,wmax,nt)
kx=np.linspace(kmin,kmax,nx)
wS1= np.sqrt(beta)*kx
wS2=-np.sqrt(beta)*kx
wR=( kx**2+np.sqrt(kx**4+4*kx**2))/2
wL=(-kx**2+np.sqrt(kx**4+4*kx**2))/2
fftrho=sf.fftshift(sf.fft2(hst[0,:,::-1]-1))
fftbz =sf.fftshift(sf.fft2(hst[5,:,::-1]))
plt.figure(figsize=(12,6))
plt.subplot(121);plt.imshow(abs(fftrho),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),cmap='Blues')
plt.plot(kx,wS1,'r--')
plt.plot(kx,wS2,'r--')
plt.axis([kmin/4,kmax/4,0,wmax/2])
plt.title(r'$\rho$');plt.xlabel(r'$k(c/\omega_{pi})$');plt.ylabel(r'$\omega/\Omega_{ci}$')
plt.clim(0,1e1)
plt.subplot(122);plt.imshow(abs(fftbz),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),cmap='Blues')
plt.plot(kx,wR,'r--')
plt.plot(kx,wL,'r--')
plt.title(r'$b_z$');plt.xlabel(r'$k(c/\omega_{pi})$');plt.ylabel(r'$\omega/\Omega_{ci}$')
plt.axis([kmin/4,kmax/4,0,wmax/2])
plt.clim(0,1e1)
plt.show()
Modulational Instability¶
We put large amplitude right-handed circularly polarized waves and white noises as an initial condition.
In [6]:
%%time
nx=1024;nt=1024*30;dx=0.25;dt=0.1*dx;isav=15
U=np.zeros((6,nx))
x=np.arange(nx)*dx
U[0,:]=1
U[4,:]=U[4,:]+0.4*np.sin(2*np.pi/(nx*dx)*x*20)
U[5,:]=U[5,:]-0.4*np.cos(2*np.pi/(nx*dx)*x*20)
U[2,:]=-(1.275)*U[4,:]
U[3,:]=-(1.275)*U[5,:]
for ik in range(256):
ikk=ik-128
#U[0,:]=U[0,:]+1e-4*np.random.randn()*np.sin(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())
U[4,:]=U[4,:]+1e-4*np.random.randn()*np.sin(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())#U[4,jm1]/4+U[4,:]/2+U[4,jp1]/4
U[5,:]=U[5,:]+1e-4*np.random.randn()*np.cos(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())#U[4,jm1]/4+U[4,:]/2+U[4,jp1]/4
bx=1; beta=0.5
data=HallMHD(nx,nt,dx,dt,beta,bx,U,isav)
locals().update(data)
Time-Space Evolution¶
In [7]:
plt.figure(figsize=(10,5))
plt.subplot(121);plt.imshow(hst[0,:,:],origin='lower',aspect='auto', extent=([0,nx*dx,0,nt*dt] ));plt.xlabel('x');plt.ylabel('t');plt.title(r'$\rho$')
plt.colorbar()
plt.subplot(122);plt.imshow(np.sqrt(hst[4,:,:]**2+hst[5,:,:]**2),origin='lower',aspect='auto',cmap='terrain',extent=([0,nx*dx,0,nt*dt]));plt.xlabel('x');plt.ylabel('t');plt.title(r'$\sqrt{b_y^2+b_z^2}$')
plt.colorbar()
plt.show()
Dispersion Relation¶
In [8]:
ddt=dt*isav
ddx=dx
nd=nt//isav
kmin=2*np.pi/(ddx*nx)*(-nx/2); kmax=2*np.pi/(ddx*nx)*(nx/2)
wmin=2*np.pi/(ddt*nd)*(-nd/2); wmax=2*np.pi/(ddt*nd)*(nd/2)
kaxis=np.linspace(kmin,kmax,nx); waxis=np.linspace(wmin,wmax,nt)
kx=np.linspace(kmin,kmax,nx)
wS=np.sqrt(beta)*kx
wR=( kx**2+np.sqrt(kx**4+4*kx**2))/2
wL=(-kx**2+np.sqrt(kx**4+4*kx**2))/2
fftrho=sf.fftshift(sf.fft2(hst[0,:,::-1]-1))
fftbz =sf.fftshift(sf.fft2(hst[5,:,::-1]))
plt.figure(figsize=(12,6))
plt.subplot(121);plt.imshow(abs(fftrho),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),cmap='Blues')
plt.axis([kmin/4,kmax/4,0,wmax/2])
plt.title(r'$\rho$')
plt.xlabel(r'$k(c/\omega_{pi})$');plt.ylabel(r'$\omega/\Omega_{ci}$')
plt.clim(0,1e4)
plt.subplot(122);plt.imshow(abs(fftbz),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),cmap='Blues')
plt.xlabel(r'$k(c/\omega_{pi})$');plt.ylabel(r'$\omega/\Omega_{ci}$')
plt.axis([kmin/4,kmax/4,0,wmax/2])
plt.title(r'$b_z$')
plt.clim(0,1e3)
plt.show()
Power Spectrum¶
In [9]:
it=800
rho =hst[0,it,:]
by =hst[4,it,:]
bz =hst[5,it,:]
dx=2*np.pi/nx;
kmin=(-nx/2); kmax=(nx/2)
kx=np.linspace(kmin,kmax,nx)
rho2fft=np.abs(sf.fftshift(sf.fft(rho)))
brfft=np.abs(sf.fftshift(sf.fft((by+1j*bz))))
plt.plot(kx,rho2fft,label=(r'$\rho$'))
plt.plot(kx,brfft,label=(r'$|B_y+iB_z|$'))
plt.loglog()
plt.legend()
plt.xlabel('mode number'); plt.ylabel('power spectrum')
plt.xlim(1,nx/3)
plt.ylim(1e-4,1e6); plt.yticks(np.logspace(-4,6,6))
plt.title('it=%.2f' %(it*dt*isav))
plt.annotate('pump wave', xy=(2.2e1, 4e2), xytext=(1e1, 1e4),arrowprops=dict(facecolor='black', shrink=0.15),horizontalalignment='right',verticalalignment='bottom')
plt.show()
In [10]:
%%writefile temp.tplx
((*- extends 'article.tplx' -*))
((* block author *))
\author{Masaru Nakanotani}
((* endblock author *))
((* block title *))
\title{Hall-MHD: Linear Waves and Modulational Instability}
((* endblock title *))
writing temp.tplx
In [15]:
!jupyter nbconvert --no-input Hall-MHD.ipynb --to pdf --template temp