Hasegawa-Wakatani Equation

reference: Wakatani1984PoF, Hasegawa1987PRL, Numata2007PoP

Basic Equations

\begin{gathered} {\partial _t}\zeta + \left\{ {\phi ,\zeta } \right\} = \alpha \left( {\phi - n} \right) - \mu {\nabla ^4}\zeta \hfill \\ {\partial _t}n + \left\{ {\phi ,n} \right\} = \alpha \left( {\phi - n} \right) - \kappa {\partial _y}\phi - \mu {\nabla ^4}n \hfill \\ \zeta \equiv {\nabla ^2}\phi \hfill \\ \end{gathered}$$\left\{ {a,b} \right\} \equiv {\partial _x}a{\partial _y}b - {\partial _y}a{\partial _x}b$$

Python Packages

%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
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['animation.embed_limit']=60
plt.rcParams['animation.html'] = 'jshtml'

Source Code

def HW(nx,ny,lx,ly,nt,dt,kap,alph,mu,phi,n,isav):
    global KX,KY,KX2,KY2,KXD,KYD
    
    dx=lx/nx; dy=ly/ny

    ### define grids ###
    kx =2*np.pi/lx*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
    ky =2*np.pi/ly*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)

    phif=sf.fft2(phi)
    nf  =sf.fft2(n)
    zetaf=-(KX2+KY2)*phif
    
    phihst =np.zeros((nt//isav,nx,ny))
    nhst   =np.zeros((nt//isav,nx,ny))
    zetahst=np.zeros((nt//isav,nx,ny))
    phihst[0,:,:] =np.real(sf.ifft2(phif))
    nhst[0,:,:]   =np.real(sf.ifft2(nf))
    zetahst[0,:,:]=np.real(sf.ifft2(zetaf))

    for it in range(1,nt):

        #---double steps with integrating factor method(4th-order Runge-Kutta)---#
        zetaf=np.exp(-mu*(KX2+KY2)**2*dt)*zetaf
        nf   =np.exp(-mu*(KX2+KY2)**2*dt)*nf

        gw1,ga1=adv(zetaf           ,nf           )
        gw2,ga2=adv(zetaf+0.5*dt*gw1,nf+0.5*dt*ga1)
        gw3,ga3=adv(zetaf+0.5*dt*gw2,nf+0.5*dt*ga2)
        gw4,ga4=adv(zetaf+    dt*gw3,nf+    dt*ga3)
        
        zetaf=zetaf+dt*(gw1+2*gw2+2*gw3+gw4)/6
        nf   =nf   +dt*(ga1+2*ga2+2*ga3+ga4)/6

        if(it%isav==0):
            phif=zetaf/(-(KX2+KY2+1e-10)); phif[0,0]=0
            phi=np.real(sf.ifft2(phif))
            n   =np.real(sf.ifft2(nf))
            zeta=np.real(sf.ifft2(zetaf))
            phihst[it//isav,:,:]=phi  
            nhst[it//isav,:,:]=n
            zetahst[it//isav,:,:]=zeta
    
    return locals()

def adv(zetaf,nf):
    phif=zetaf/(-(KX2+KY2+1e-10)); phif[0,0]=0

    phi=np.real(sf.ifft2(phif))
    n  =np.real(sf.ifft2(nf))

    phixf = 1j*KX*phif;  phix =np.real(sf.ifft2(phixf *KXD*KYD))
    phiyf = 1j*KY*phif;  phiy =np.real(sf.ifft2(phiyf *KXD*KYD))
    zetaxf= 1j*KX*zetaf; zetax=np.real(sf.ifft2(zetaxf*KXD*KYD))
    zetayf= 1j*KY*zetaf; zetay=np.real(sf.ifft2(zetayf*KXD*KYD))
    nxf   = 1j*KX*nf;    nnx   =np.real(sf.ifft2(nxf   *KXD*KYD))
    nyf   = 1j*KY*nf;    nny   =np.real(sf.ifft2(nyf   *KXD*KYD))
    
    advf =-(phix*zetay-phiy*zetax)+alph*(phi-n)
    advg =-(phix*nny  -phiy*nnx)  +alph*(phi-n)-kap*np.real(sf.ifft2(phiyf))
    advff=sf.fft2(advf)
    advgf=sf.fft2(advg)
    
    return advff,advgf

Initial Condition

nx=256; ny=256; nt=500; isav=50
kap=1.0
alph=0.1
mu=1e-4
dt=2e-2
lx=2*np.pi/0.15; ly=2*np.pi/0.15
dx=lx/nx; dy=ly/ny
x  =np.arange(nx)*dx
y  =np.arange(ny)*dy
X,Y=np.meshgrid(x,y)

s=2; s2=s**2
r1=(X-lx/2)**2+(Y-ly/2)**2
n  =np.exp(-r1/s2)
phi=n

Run

%%time
data=HW(nx,ny,lx,ly,nt,dt,kap,alph,mu,phi,n,isav)
locals().update(data)

CPU times: user 20.1 s, sys: 119 ms, total: 20.2 s
Wall time: 20.5 s

Animation

def update_anim(it):
    
    fig.clf()

    ax1 = fig.add_subplot(221)
    ax2 = fig.add_subplot(222)
    ax3 = fig.add_subplot(223)
    ax4 = fig.add_subplot(224)    
    
    for ax in (ax1, ax2, ax3, ax4):
        ax.clear()   

    im1=ax1.imshow(zetahst[it,:,:]            ,aspect='auto',origin='lower',cmap='jet');ax1.axis('off');fig.colorbar(im1, ax=ax1);ax1.set_title(r'$\zeta\ (vorticity)$')
    im2=ax2.imshow(nhst[it,:,:]               ,aspect='auto',origin='lower',cmap='jet');ax2.axis('off');fig.colorbar(im2, ax=ax2);ax2.set_title(r'$n\ (density)$')
    im3=ax3.imshow(phihst[it,:,:]             ,aspect='auto',origin='lower',cmap='jet');ax3.axis('off');fig.colorbar(im3, ax=ax3);ax3.set_title(r'$\phi\ (potential)$')
    im4=ax4.imshow(phihst[it,:,:]-nhst[it,:,:],aspect='auto',origin='lower',cmap='jet');ax4.axis('off');fig.colorbar(im4, ax=ax4);ax4.set_title(r'$\phi-n$')

fig=plt.figure(figsize=(10,8))
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)    
plt.close()
anim

Once Loop Reflect

Modified Hasegawa-Wakatani Equation

\begin{gathered} {\partial _t}\zeta + \left\{ {\phi ,\zeta } \right\} = \alpha \left( {\tilde \phi - \tilde n} \right) - \mu {\nabla ^4}\zeta \hfill \\ {\partial _t}n + \left\{ {\phi ,n} \right\} = \alpha \left( {\tilde \phi - \tilde n} \right) - \kappa {\partial _y}\phi - \mu {\nabla ^4}n \hfill \\ \zeta \equiv {\nabla ^2}\phi \hfill \\ \left\{ {a,b} \right\} \equiv {\partial _x}a{\partial _y}b - {\partial _y}a{\partial _x}b \hfill \\ \tilde f \equiv f - \left\langle f \right\rangle \hfill \\ \left\langle f \right\rangle \equiv \frac{1}{{{L_y}}}\int {fdy} \hfill \\ \end{gathered}

Source Code

def MHW(nx,ny,lx,ly,nt,dt,kap,alph,mu,phi,n,isav):
    global KX,KY,KX2,KY2,KXD,KYD
    
    dx=lx/nx; dy=ly/ny

    ### define grids ###
    kx =2*np.pi/lx*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
    ky =2*np.pi/ly*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)

    phif=sf.fft2(phi)
    nf  =sf.fft2(n)
    #zetaf=sf.fft2(np.random.randn(nx,ny)*2)
    zetaf=-(KX2+KY2)*phif
    
    phihst =np.zeros((nt//isav,nx,ny))
    nhst   =np.zeros((nt//isav,nx,ny))
    zetahst=np.zeros((nt//isav,nx,ny))
    phihst[0,:,:] =np.real(sf.ifft2(phif))
    nhst[0,:,:]   =np.real(sf.ifft2(nf))
    zetahst[0,:,:]=np.real(sf.ifft2(zetaf))

    for it in range(1,nt):

        #---double steps with integrating factor method(4th-order Runge-Kutta)---#
        zetaf=np.exp(-mu*(KX2+KY2)**2*dt)*zetaf
        nf   =np.exp(-mu*(KX2+KY2)**2*dt)*nf

        gw1,ga1=adv(zetaf           ,nf           )
        gw2,ga2=adv(zetaf+0.5*dt*gw1,nf+0.5*dt*ga1)
        gw3,ga3=adv(zetaf+0.5*dt*gw2,nf+0.5*dt*ga2)
        gw4,ga4=adv(zetaf+    dt*gw3,nf+    dt*ga3)
        
        zetaf=zetaf+dt*(gw1+2*gw2+2*gw3+gw4)/6
        nf   =nf   +dt*(ga1+2*ga2+2*ga3+ga4)/6

        if(it%isav==0):
            phif=zetaf/(-(KX2+KY2+1e-10)); phif[0,0]=0
            phi=np.real(sf.ifft2(phif))
            n   =np.real(sf.ifft2(nf))
            zeta=np.real(sf.ifft2(zetaf))
            phihst[it//isav,:,:]=phi  
            nhst[it//isav,:,:]=n
            zetahst[it//isav,:,:]=zeta
    
    return locals()

def adv(zetaf,nf):
    phif=zetaf/(-(KX2+KY2+1e-10)); phif[0,0]=0

    phi=np.real(sf.ifft2(phif))
    n  =np.real(sf.ifft2(nf))
    
    phiz=np.sum(phi*dy,axis=0)/ly
    nz  =np.sum(n  *dy,axis=0)/ly

    phixf = 1j*KX*phif;  phix =np.real(sf.ifft2(phixf *KXD*KYD))
    phiyf = 1j*KY*phif;  phiy =np.real(sf.ifft2(phiyf *KXD*KYD))
    zetaxf= 1j*KX*zetaf; zetax=np.real(sf.ifft2(zetaxf*KXD*KYD))
    zetayf= 1j*KY*zetaf; zetay=np.real(sf.ifft2(zetayf*KXD*KYD))
    nxf   = 1j*KX*nf;    nnx   =np.real(sf.ifft2(nxf  *KXD*KYD))
    nyf   = 1j*KY*nf;    nny   =np.real(sf.ifft2(nyf  *KXD*KYD))
    
    advf =-(phix*zetay-phiy*zetax)+alph*((phi-phiz)-(n-nz))
    advg =-(phix*nny  -phiy*nnx)  +alph*((phi-phiz)-(n-nz))-kap*np.real(sf.ifft2(phiyf))
    advff=sf.fft2(advf)
    advgf=sf.fft2(advg)
    
    return advff,advgf

Initial Condition

nx=256; ny=256; nt=5000; isav=500
kap=1.0
alph=1.0
mu=1e-4
dt=2e-2
lx=2*np.pi/0.15; ly=2*np.pi/0.15
dx=lx/nx; dy=ly/ny
x  =np.arange(nx)*dx
y  =np.arange(ny)*dy
X,Y=np.meshgrid(x,y)

s=2; s2=s**2
r1=(X-lx/2)**2+(Y-ly/2)**2
n  =np.exp(-r1/s2)
phi=n

Run

%%time
data=MHW(nx,ny,lx,ly,nt,dt,kap,alph,mu,phi,n,isav)
locals().update(data)

CPU times: user 3min 27s, sys: 482 ms, total: 3min 28s
Wall time: 3min 30s

Animation

fig=plt.figure(figsize=(10,8))
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)    
plt.close()
anim

Once Loop Reflect

!jupyter nbconvert --to html --template tmp.tpl Hasegawa-Wakatani_equation.ipynb

[NbConvertApp] Converting notebook Hasegawa-Wakatani_equation.ipynb to html
[NbConvertApp] Writing 6173892 bytes to Hasegawa-Wakatani_equation.html