Reference: Biskamp1997PoP, Biskamp1995PRL, Biskamp1999PoP
$${\partial _t}\left( {1 - d_e^2{\nabla ^2}} \right){\bf{B}} - \nabla \times \left[ {{{\bf{v}}_e} \times \left( {1 - d_e^2{\nabla ^2}} \right){\bf{B}}} \right] = -{\eta _\nu }{\left( { - {\nabla ^2}} \right)^\nu }{\bf{B}}$$\begin{array}{l} {\bf{B}} = {{\bf{e}}_z} \times \nabla \psi + {{\bf{e}}_z}\left( {{B_{z0}} + b} \right),\\ {\bf{j}} = \nabla b \times {{\bf{e}}_z} + {{\bf{e}}_z}{\nabla ^2}\psi = - {{\bf{v}}_e} \end{array}\begin{array}{l} {\partial _t}\left( {\psi - d_e^2j_z} \right) + {{\bf{v}}_e} \cdot \nabla \left( {\psi - d_e^2j_z} \right) = -{\eta _\nu }{\left( { - {\nabla ^2}} \right)^\nu }\psi ,\\ {\partial _t}\left( {b - d_e^2w} \right) + {{\bf{v}}_e} \cdot \nabla \left( {b - d_e^2w} \right) + {\bf{B}} \cdot \nabla j_z = -{\eta _\nu }{\left( { - {\nabla ^2}} \right)^\nu }b,\\ j_z = {\nabla ^2}\psi ,\\ w = {\nabla ^2}b \end{array}%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 numba.decorators import jit
from tqdm import tqdm
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
def EMHD(nx,ny,nt,dt,de2,eta,nu,psif,bf,isav):
global KX,KY,KX2,KY2,KXD,KYD
dx=2*np.pi/nx; dy=2*np.pi/ny
### define grids ###
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
KX ,KY =np.meshgrid(kx ,ky )
KX2,KY2=np.meshgrid(kx2,ky2)
KXD,KYD=np.meshgrid(kxd,kyd)
bhst =np.zeros((nt//isav,nx,ny))
psihst=np.zeros((nt//isav,nx,ny))
jzhst =np.zeros((nt//isav,nx,ny))
whst =np.zeros((nt//isav,nx,ny))
jz =np.real(sf.ifft2(-(KX2+KY2)*psif))
w =np.real(sf.ifft2(-(KX2+KY2)*bf))
Ehst =np.zeros(nt)
psihst[0,:,:]=np.real(sf.ifft2(psif))
jzhst[0,:,:]=jz
bhst[0,:,:]=np.real(sf.ifft2(bf))
whst[0,:,:]=w
Ehst[0]=np.sum(np.real(sf.ifft2(1j*KX*psif))**2+np.real(sf.ifft2(1j*KY*psif))**2+np.real(sf.ifft2(bf))**2+de2*(np.real(sf.ifft2(-(KX**2+KY**2)*psif))**2+np.real(sf.ifft2(1j*KX*bf))**2+np.real(sf.ifft2(1j*KY*bf))**2))
for it in range(1,nt):
#---double steps with integrating factor method(4th-order Runge-Kutta)---#
psif=np.exp(-eta*(KX2+KY2)**nu/(1+de2*(KX2+KY2))*dt)*psif
bf =np.exp(-eta*(KX2+KY2)**nu/(1+de2*(KX2+KY2))*dt)*bf
jzf =-(KX2+KY2)*psif
wf =-(KX2+KY2)*bf
ff =psif-de2*jzf
gf =bf -de2*wf
#---4th-order Runge-Kutta---#
gw1,ga1=adv(ff ,gf )
gw2,ga2=adv(ff+0.5*dt*gw1,gf+0.5*dt*ga1)
gw3,ga3=adv(ff+0.5*dt*gw2,gf+0.5*dt*ga2)
gw4,ga4=adv(ff+ dt*gw3,gf+ dt*ga3)
ff=ff+dt*(gw1+2*gw2+2*gw3+gw4)/6.0
gf=gf+dt*(ga1+2*ga2+2*ga3+ga4)/6.0
psif=ff/(1.0+de2*(KX2+KY2))
bf =gf/(1.0+de2*(KX2+KY2))
Ehst[it]=np.sum(np.real(sf.ifft2(KX*psif))**2+np.real(sf.ifft2(KY*psif))**2+np.real(sf.ifft2(bf))**2+de2*(np.real(sf.ifft2((KX**2+KY**2)*psif))**2+np.real(sf.ifft2(KX*bf))**2+np.real(sf.ifft2(KY*bf))**2))
if(it%isav==0):
psi=np.real(sf.ifft2(psif))
b =np.real(sf.ifft2(bf))
jz =np.real(sf.ifft2(-(KX2+KY2)*psif))
w =np.real(sf.ifft2(-(KX2+KY2)*bf))
psihst[it//isav,:,:]=psi
jzhst[it//isav,:,:]=jz
bhst[it//isav,:,:]=b
whst[it//isav,:,:]=w
#plt.figure(figsize=(8,8))
#kxmin=2*np.pi/(dx*nx)*(-nx/2); kxmax=2*np.pi/(dx*nx)*(nx/2)
#kymin=2*np.pi/(dy*ny)*(-ny/2); kymax=2*np.pi/(dy*ny)*(ny/2)
#Kp=np.sqrt(KX**2+KY**2)
#fftv2 =np.abs(bf)**2+np.abs(-1j*KY*psif)**2+np.abs(1j*KX*psif)**2
#(spc,bns) =np.histogram(Kp,bins=nx,weights=fftv2)
#binx=0.5*(bns[1:]+bns[:-1])
#display.clear_output(True)
#plt.plot(binx,spc,'.-')
#plt.ylim(1e-10,1e11)
##plt.plot(kx,kx**(7/3)*fftv2[0,:])
#plt.title(r'$it=%05d$'%(it))
#plt.xlim(1,1e3)
#plt.xlabel('k')
#plt.ylabel(r'$P(k)$')
#plt.loglog()
#plt.show()
return locals()
def adv(ff,gf):
ff=ff*KXD*KYD
gf=gf*KXD*KYD
psif=ff/(1+de2*(KX2+KY2))
jzf =-(KX2+KY2)*psif
bf =gf/(1+de2*(KX2+KY2))
vexf=-1j*KY*bf ; vex=np.real(sf.ifft2(vexf))
veyf= 1j*KX*bf ; vey=np.real(sf.ifft2(veyf))
Bxf =-1j*KY*psif; Bx =np.real(sf.ifft2(Bxf ))
Byf = 1j*KX*psif; By =np.real(sf.ifft2(Byf ))
fxf = 1j*KX*ff ; fx =np.real(sf.ifft2(fxf ))
fyf = 1j*KY*ff ; fy =np.real(sf.ifft2(fyf ))
gxf = 1j*KX*gf ; gx =np.real(sf.ifft2(gxf ))
gyf = 1j*KY*gf ; gy =np.real(sf.ifft2(gyf ))
jzxf= 1j*KX*jzf ; jzx=np.real(sf.ifft2(jzxf))
jzyf= 1j*KY*jzf ; jzy=np.real(sf.ifft2(jzyf))
advf =vex*fx+vey*fy
advg =vex*gx+vey*gy
advj =Bx*jzx+By*jzy
advff=sf.fft2(advf)
advgf=sf.fft2(advg)
advjf=sf.fft2(advj)
return -advff,-advgf-advjf
%%time
nx=512; ny=512; nt=2000; isav=10
nu=3
de=1; de2=de**2
eta=4e-8
dt=2e-3
dx=2*np.pi/nx; dy=2*np.pi/ny
x =np.arange(nx)*dx
y =np.arange(ny)*dy
X,Y=np.meshgrid(x,y)
psi=np.cos(2*X+2.3)+np.cos(Y+4.1)
b =np.cos( X+1.4)+np.cos(Y+0.5)
psif=sf.fft2(psi)
bf =sf.fft2(b)
data=EMHD(nx,ny,nt,dt,de2,eta,nu,psif,bf,isav)
locals().update(data)
plt.rcParams['animation.html'] = 'html5'
fig=plt.figure(figsize=(10,10))
ax1 = fig.add_subplot(111)
def update_anim(it):
ax1.clear()#; ax2.clear()
ax1.imshow(whst[it,:,:],origin='lower',aspect='auto',cmap='gray')
ax1.axis('off')
ax1.set_title(r'$w(x,y)=\nabla^2b$ (t=%.2f)' %(it*dt*isav))
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)
plt.close()
anim