MHD Simulation

This post shows MHD simulations including linear waves and the parametric decay instability.

Governing Equations

We consider 1D case and assume the equation of state to be isothermal:

\begin{equation} \begin{gathered} {\partial _t}\rho + {\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} - vb_x} \right); \hfill \\ {\partial _t}{b_z} = - {\partial _x}\left( {u{b_z} - wb_x} \right); \hfill \\ p = \beta \rho . \hfill \\ \end{gathered} \end{equation}
In [3]:
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as sf
import matplotlib.animation as animation
plt.rcParams['font.size'] = 12
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['animation.html'] = 'jshtml'
In [6]:
def MHD(nx,nt,dx,dt,beta,bx,U,isav):
    global kx,kd
    
    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)

    for it in range(1,nt):
        
        #---add viscocity and resistivity term---#
        mu =1e-5*(1+(kx/4)**2)
        eta=1e-5*(1+(kx/4)**2)
        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')
    ro =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(ro*u)
    ff[1,:]=-sf.fft(u*ux+by*byx/ro+bz*bzx/ro)-1j*kx*sf.fft(np.log(ro))*beta
    ff[2,:]=-sf.fft(u*vx-bx*byx/ro)
    ff[3,:]=-sf.fft(u*wx-bx*bzx/ro)
    ff[4,:]=-1j*kx*sf.fft(u*by-v*bx)
    ff[5,:]=-1j*kx*sf.fft(u*bz-w*bx)
 
    return ff

Linear Waves

Linearizing and Fourier transforming, we obtain the dispersion realtion of the system.
These are Alfven ($\omega=V_A k$) and acoustic modes ($\omega=C_Sk$).

In [7]:
%%time
nx=2048;nt=512*40;dx=2*np.pi/nx;dt=0.1*dx;isav=40
U=np.zeros((6,nx))
x=np.arange(nx)*dx

U[0,:]=1

for ik in range(128):
    ikk=ik-64
    U[0,:]=U[0,:]+1e-5*np.random.randn()*np.sin(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())
    U[4,:]=U[4,:]+1e-5*np.random.randn()*np.cos(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())
    U[5,:]=U[5,:]+1e-5*np.random.randn()*np.cos(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())

bx=1; beta=0.1
data=MHD(nx,nt,dx,dt,beta,bx,U,isav)
locals().update(data)

CPU times: user 44.1 s, sys: 93.6 ms, total: 44.1 s
Wall time: 44.6 s

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)
wS1= np.sqrt(beta)*kx
wS2=-np.sqrt(beta)*kx
wR= kx
wL=-kx

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/16,kmax/16,0,wmax/2])
plt.title(r'$\rho$')
plt.xlabel(r'$k$')
plt.ylabel(r'$\omega$')
plt.clim(0,1e0)

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$')
plt.ylabel(r'$\omega$')
plt.axis([kmin/16,kmax/16,0,wmax/2])
plt.clim(0,1e0)
plt.show()

Parametric Instability of a Large Amplitude Alfven Wave

We set a large amplitude left-handed circularly polarized wave and white noises in the background plasma as an initial condition. The Alfven mode experiences the parametric instability and is decayed into a backward Alfven wave and an acoustic wave matching the conditions, $k_1=k_2+k_3$ and $\omega_1=\omega_2+\omega_3$.

In [13]:
%%time
nx=2048;nt=512*250;dx=2*np.pi/nx;dt=0.1*dx;isav=250
U=np.zeros((6,nx))
x=np.arange(nx)*dx

U[0,:]=1
U[4,:]= 0.2*np.cos(4*x)
U[5,:]=-0.2*np.sin(4*x)
U[2,:]=-(1.0)*U[4,:]
U[3,:]=-(1.0)*U[5,:]

for ik in range(128):
    ikk=ik-64
    U[0,:]=U[0,:]+1e-5*np.random.randn()*np.sin(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())
    U[4,:]=U[4,:]+1e-5*np.random.randn()*np.sin(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())
    U[5,:]=U[5,:]+1e-5*np.random.randn()*np.cos(2*np.pi/(nx*dx)*x*ikk+2*np.pi*np.random.rand())

bx=1; beta=0.1
data=MHD(nx,nt,dx,dt,beta,bx,U,isav)
locals().update(data)

CPU times: user 4min 33s, sys: 448 ms, total: 4min 34s
Wall time: 4min 38s
In [14]:
### Fourier decomposition (Terasawa1986JGR) 
lx=nx*dx
x=np.linspace(-lx/2,lx/2,nx)
byc=np.zeros((nt//isav,nx//2))
bys=np.zeros((nt//isav,nx//2))
bzc=np.zeros((nt//isav,nx//2))
bzs=np.zeros((nt//isav,nx//2))
for it in range(nt//isav):
    for ik in range(nx//2):
        byc[it,ik]=2/lx*np.sum(hst[4,it,:]*np.cos(2*np.pi*(ik+1)*x/lx)*dx)
        bys[it,ik]=2/lx*np.sum(hst[4,it,:]*np.sin(2*np.pi*(ik+1)*x/lx)*dx)
        bzc[it,ik]=2/lx*np.sum(hst[5,it,:]*np.cos(2*np.pi*(ik+1)*x/lx)*dx)
        bzs[it,ik]=2/lx*np.sum(hst[5,it,:]*np.sin(2*np.pi*(ik+1)*x/lx)*dx)

fftbr=0.5*((byc+bzs)+1j*(bzc-bys))
fftbl=0.5*((byc-bzs)+1j*(bzc+bys))

br=np.zeros((nt//isav,nx),dtype='complex')
bl=np.zeros((nt//isav,nx),dtype='complex')
for it in range(nt//isav):
    for ik in range(nx//2):
        br[it,:]=br[it,:]+fftbr[it,ik]*np.exp( 1j*2*np.pi*(ik+1)*x/lx)
        bl[it,:]=bl[it,:]+fftbl[it,ik]*np.exp(-1j*2*np.pi*(ik+1)*x/lx)

byr=np.real(br)
bzr=np.imag(br)
byl=np.real(bl)
bzl=np.imag(bl)
In [15]:
def update_anim(it):
    
    fig.clf()

    ax1 = fig.add_subplot(111)
    
    ax1.clear() 

    xx=dx*np.arange(nx)
    im1=ax1.plot(xx,hst[0,it*4,:]-1,label=r'$\rho$')
    im2=ax1.plot(xx,byr[it*4,:],label=r'$b_y\;(\leftarrow L^{-}\; \;R^{+}\rightarrow)$')
    im3=ax1.plot(xx,byl[it*4,:],label=r'$b_y\;(\leftarrow R^{-}\; \;L^{+}\rightarrow)$')
    ax1.legend(loc='upper right')
    ax1.set_xlabel(r'$X$')
    ax1.set_ylabel(r'$\Psi$')
    ax1.set_xlim(0,nx*dx)
    ax1.set_ylim(-1,1)
In [16]:
fig=plt.figure(figsize=(16,5))
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav//4)
plt.close()
anim
Out[16]:
In [18]:
xx=dx*np.arange(nx)
fig, ax = plt.subplots(ncols=3, sharex=True, figsize=(16, 8))
ax[0].imshow(hst[0,:,:],cmap='jet',origin='lower',aspect='auto',extent=([0,nx*dx,0,nt*dt]))
ax[1].imshow(byl[:,:]  ,cmap='bwr',origin='lower',aspect='auto',extent=([0,nx*dx,0,nt*dt]),vmin=-0.2,vmax=0.2)
ax[2].imshow(byr[:,:]  ,cmap='bwr',origin='lower',aspect='auto',extent=([0,nx*dx,0,nt*dt]),vmin=-0.2,vmax=0.2)
ax[0].set_title(r'$\rho$')
ax[1].set_title(r'$b_y\;(\leftarrow R^{-}\; \;L^{+}\rightarrow)$')
ax[2].set_title(r'$b_y\;(\leftarrow L^{-}\; \;R^{+}\rightarrow)$')
for i in range(3):
    ax[i].set_xlabel(r'$X$')
ax[0].set_ylabel(r'$T$')
plt.show()

Dispersion Relation

In [19]:
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
wL=-kx

fftrho=sf.fftshift(sf.fft2(hst[0,:,::-1]-1))
fftbz =sf.fftshift(sf.fft2(hst[5,:,::-1]))

plt.figure(figsize=(12,6))

plt.imshow(abs(fftrho)+abs(fftbz),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),cmap='Blues')
plt.plot(kx,wR ,'k--',lw=1)
plt.plot(kx,wL ,'k--',lw=1)
plt.plot(kx,wS1,'k--',lw=1)
plt.plot(kx,wS2,'k--',lw=1)
plt.axis([kmin/32,kmax/32,0,wmax/4])
plt.xlabel(r'$k$')
plt.ylabel(r'$\omega$')
plt.clim(0,2e4)

plt.arrow(0,0,4,4, head_width=0.01, head_length=0.01, fc='C1', ec='C1',zorder=5)
plt.arrow(0,0,6,np.sqrt(beta)*6,head_width=0.01, head_length=0.01, fc='C1', ec='C1',zorder=5)
plt.arrow(0,0,-2,2,head_width=0.01, head_length=0.01, fc='C1', ec='C1',zorder=5)
plt.arrow(4,4,-6,-np.sqrt(beta)*6,head_width=0.01, head_length=0.01, fc='C1', ec='C1',ls='--',width=0.001)
plt.arrow(4,4,2,-2,head_width=0.01, head_length=0.01, fc='C1', ec='C1',ls='--',width=0.001)
plt.show()

Higher harmonics on the ascoutic mode branch can be understood as the result of the wave steepening.

In [14]:
!jupyter nbconvert --to html --template tmp.tpl MHD_parametric_instability.ipynb

[NbConvertApp] Converting notebook MHD_parametric_instability.ipynb to html
[NbConvertApp] Writing 6690929 bytes to MHD_parametric_instability.html