\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}
%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
from numba.decorators import jit
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
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.r_[np.ones(nx//3),np.zeros(nx//3+nx%3),np.ones(nx//3)] #de-aliased
kyd=np.r_[np.ones(ny//3),np.zeros(ny//3+ny%3),np.ones(ny//3)] #de-aliased
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);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);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
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)
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,:,:],80,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