Reference: 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}
\begin{array}{l} \psi \left( {x,y} \right) = \cos \left( {2x + 2.3} \right) + \cos \left( {y + 4.1} \right)\\ b\left( {x,y} \right) = \cos \left( {x + 1.4} \right) + \cos \left( {y + 0.5} \right) \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,psi,b,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.r_[np.ones(nx//3),np.zeros(nx//3+nx%3),np.ones(nx//3)] #for de-aliasing
kyd=np.r_[np.ones(ny//3),np.zeros(ny//3+ny%3),np.ones(ny//3)] #for de-aliasing
kx2=kx**2; ky2=ky**2
KX ,KY =np.meshgrid(kx ,ky )
KX2,KY2=np.meshgrid(kx2,ky2)
KXD,KYD=np.meshgrid(kxd,kyd)
psif=sf.fft2(psi)
bf =sf.fft2(b)
bhst =np.zeros((nt//isav,nx,ny))
psihst=np.zeros((nt//isav,nx,ny))
whst =np.zeros((nt//isav,nx,ny))
w =np.real(sf.ifft2(-(KX2+KY2)*bf))
psihst[0,:,:]=psi
bhst[0,:,:]=b
whst[0,:,:]=w
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
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
gf=gf+dt*(ga1+2*ga2+2*ga3+ga4)/6
psif=ff/(1+de2*(KX2+KY2))
bf =gf/(1+de2*(KX2+KY2))
if(it%isav==0):
psi=np.real(sf.ifft2(psif))
b =np.real(sf.ifft2(bf))
w =np.real(sf.ifft2(-(KX2+KY2)*bf))
psihst[it//isav,:,:]=psi
bhst[it//isav,:,:]=b
whst[it//isav,:,:]=w
return locals()
def adv(ff,gf):
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*KXD*KYD))
veyf= 1j*KX*bf ; vey=np.real(sf.ifft2(veyf*KXD*KYD))
Bxf =-1j*KY*psif; Bx =np.real(sf.ifft2(Bxf *KXD*KYD))
Byf = 1j*KX*psif; By =np.real(sf.ifft2(Byf *KXD*KYD))
fxf = 1j*KX*ff ; fx =np.real(sf.ifft2(fxf *KXD*KYD))
fyf = 1j*KY*ff ; fy =np.real(sf.ifft2(fyf *KXD*KYD))
gxf = 1j*KX*gf ; gx =np.real(sf.ifft2(gxf *KXD*KYD))
gyf = 1j*KY*gf ; gy =np.real(sf.ifft2(gyf *KXD*KYD))
jzxf= 1j*KX*jzf ; jzx=np.real(sf.ifft2(jzxf*KXD*KYD))
jzyf= 1j*KY*jzf ; jzy=np.real(sf.ifft2(jzyf*KXD*KYD))
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=1e-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)
data=EMHD(nx,ny,nt,dt,de2,eta,nu,psi,b,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
I must note that the picture at $t=2.0$ above is different with Fig. 1(a) in Biskamp1999PoP. Unfortunatelly, I don't see the reason. The difference, however, would be not big because the entire devcelopment above is quite similiar.