Vorticity Transport Equation in RMHD

\begin{gathered} \frac{{\partial \omega }}{{\partial t}} = - \left( {{\mathbf{v}} \cdot \nabla } \right)\omega + \left( {{\mathbf{b}} \cdot \nabla } \right)j + R_\nu ^{ - 1}{\nabla ^2}\omega \hfill \\ \frac{{\partial a}}{{\partial t}} = - \left( {{\mathbf{v}} \cdot \nabla } \right)a + R_\mu ^{ - 1}{\nabla ^2}a \hfill \\ \end{gathered}
In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as sf
from IPython import display
import math as mt
import matplotlib.animation as animation
from tqdm import tqdm
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [2]:
def RMHD(nx,ny,nt,dt,mu,nu,w,j,isav):
    global dx,KX,KY,KX2,KY2,KXD,KYD
    
    dx=2*np.pi/nx; dy=2*np.pi/ny

    ### define grids ###
    x  =np.arange(nx)*dx
    y  =np.arange(ny)*dy
    kx =2*np.pi/(nx*dx)*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
    ky =2*np.pi/(ny*dy)*np.r_[np.arange(ny/2),np.arange(-ny/2,0)]
    kxd=np.zeros(nx)
    kyd=np.zeros(ny)
    kxd[np.abs(kx)<=2/3*np.max(abs(kx))]=1.0 ### 2/3-rule for dealiasing
    kyd[np.abs(ky)<=2/3*np.max(abs(ky))]=1.0 ### 2/3-rule for dealiasing
    kx2=kx**2; ky2=ky**2
    X  ,Y  =np.meshgrid(x,y)
    KX ,KY =np.meshgrid(kx ,ky )
    KX2,KY2=np.meshgrid(kx2,ky2)
    KXD,KYD=np.meshgrid(kxd,kyd)
    
    #w = np.exp(-((X-np.pi+np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.3)\
    #   -np.exp(-((X-np.pi-np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.2)\
    #   +np.exp(-((X-np.pi-np.pi/5)**2+(Y-np.pi-np.pi/5)**2)/0.4)
    #j = np.exp(-((X-np.pi+np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.3)\
    #   +np.exp(-((X-np.pi-np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.2)
    wf=sf.fft2(w)
    jf=sf.fft2(j)
    af=jf/(KX2+KY2+1e-10);af[0,0]=0.0

    whst=np.zeros((nt//isav,nx,ny))
    ahst=np.zeros((nt//isav,nx,ny))
    jhst=np.zeros((nt//isav,nx,ny))
    for it in tqdm(range(nt)):

        #---double steps with integrating factor method(4th-order Runge-Kutta)---#
        wf=np.exp(-nu*(KX2+KY2)*dt)*wf
        af=np.exp(-mu*(KX2+KY2)*dt)*af
        gw1,ga1=f(wf           ,af           )
        gw2,ga2=f(wf+0.5*dt*gw1,af+0.5*dt*ga1)
        gw3,ga3=f(wf+0.5*dt*gw2,af+0.5*dt*ga2)
        gw4,ga4=f(wf+    dt*gw3,af+    dt*ga3)
        wf=wf+dt*(gw1+2*gw2+2*gw3+gw4)/6
        af=af+dt*(ga1+2*ga2+2*ga3+ga4)/6

        if(it%isav==0):
            w=np.real(sf.ifft2(wf))
            a=np.real(sf.ifft2(af))
            j=np.real(sf.ifft2(af*(KX2+KY2)))
            whst[it//isav,:,:]=w  
            ahst[it//isav,:,:]=a
            jhst[it//isav,:,:]=j
    
    return locals()

def f(wf,af):
    jf  =af*(KX2+KY2)
    phif=wf/(KX2+KY2+1e-10);phif[0,0]=0.0
    vxf = 1j*KY*phif; vx=np.real(sf.ifft2(vxf*KXD*KYD))
    vyf =-1j*KX*phif; vy=np.real(sf.ifft2(vyf*KXD*KYD))
    bxf = 1j*KY*af  ; bx=np.real(sf.ifft2(bxf*KXD*KYD))
    byf =-1j*KX*af  ; by=np.real(sf.ifft2(byf*KXD*KYD))
    wxf = 1j*KX*wf  ; wx=np.real(sf.ifft2(wxf*KXD*KYD))
    wyf = 1j*KY*wf  ; wy=np.real(sf.ifft2(wyf*KXD*KYD))
    axf = 1j*KX*af  ; ax=np.real(sf.ifft2(axf*KXD*KYD))
    ayf = 1j*KY*af  ; ay=np.real(sf.ifft2(ayf*KXD*KYD))
    jxf = 1j*KX*jf  ; jx=np.real(sf.ifft2(jxf*KXD*KYD))
    jyf = 1j*KY*jf  ; jy=np.real(sf.ifft2(jyf*KXD*KYD))
    advw =vx*wx+vy*wy
    advj =bx*jx+by*jy
    adva =vx*ax+vy*ay
    advwf=sf.fft2(advw)
    advjf=sf.fft2(advj)
    advaf=sf.fft2(adva)

    return -advwf+advjf,-advaf
In [3]:
nx=512; ny=512; nt=1200; isav=10
mu=0.0002; nu=0.0002
dt=0.08
w=np.random.randn(nx,ny)*5
j=np.random.randn(nx,ny)*5

data=RMHD(nx,ny,nt,dt,mu,nu,w,j,isav)
locals().update(data)

100%|█████████████████████████████████| 1200/1200 [04:37<00:00,  4.32it/s]
In [4]:
plt.rcParams['animation.html'] = 'html5'
fig=plt.figure(figsize=(8,8))
ax1 = fig.add_subplot(111)

def update_anim(it):    
    ax1.clear()#; ax2.clear()

    ax1.imshow(jhst[it,:,:],origin='lower',aspect='auto',cmap='jet')
    ax1.contour(ahst[it,:,:],40,linestyles='solid',linewidths=0.5,colors='k')
    ax1.set_title(r'$j(color\ map)+a(contour\ plot)$')
    ax1.set_xlabel(r'$x$')
    ax1.set_ylabel(r'$y$')

anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)    
plt.close()
anim
Out[4]: