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

/var/folders/hl/2w9c9jdj5jqgbpmm44p6whpr0000gn/T/ipykernel_72160/1714137282.py:2: DeprecationWarning: Importing display from IPython.core.display is deprecated since IPython 7.14, please import from IPython display
  from IPython.core.display import display, HTML, clear_output
/var/folders/hl/2w9c9jdj5jqgbpmm44p6whpr0000gn/T/ipykernel_72160/1714137282.py:2: DeprecationWarning: Importing clear_output from IPython.core.display is deprecated since IPython 7.14, please import from IPython display
  from IPython.core.display import display, HTML, clear_output
In [13]:
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.zeros(nx)
    kd[np.abs(kx)<=2/3*np.max(abs(kx))]=1.0 ### 2/3-rule for dealiasing
 
    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 [14]:
%%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)

CPU times: user 57.7 s, sys: 169 ms, total: 57.8 s
Wall time: 59.3 s

Time-Space Evolution

In [15]:
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 [16]:
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,kmax,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 [17]:
%%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)

CPU times: user 44 s, sys: 179 ms, total: 44.2 s
Wall time: 45.7 s

Time-Space Evolution

In [18]:
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((hst[5,:,:]),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()
plt.plot(hst[0,:,0])
Out[18]:
[<matplotlib.lines.Line2D at 0x7f8dc0ad1b20>]

Dispersion Relation

In [20]:
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 [21]:
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 [ ]: