Vorticity Transport Equation

$$\frac{{\partial \omega }}{{\partial t}} = - \left( {{\mathbf{v}} \cdot \nabla } \right)\omega + R_\nu ^{ - 1}{\nabla ^2}\omega$$

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
from numba.decorators import jit
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [2]:
def VTE(nx,ny,nt,dt,mu,w,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
    X,Y=np.meshgrid(x,y)
    kx =2*np.pi/(dx*nx)*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
    ky =2*np.pi/(dy*ny)*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
    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/10)**2+(Y-np.pi)**2)/0.1)\
    #   +np.exp(-((X-np.pi-np.pi/10)**2+(Y-np.pi)**2)/0.1)
    #   +np.exp(-((X-2*np.pi-np.pi/10)**2+(Y-2*np.pi-np.pi/5)**2)/0.1)
    wf=sf.fft2(w)

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

        ### Cranck-Nicholson (only for the dissipation term)###
        #wf=1/(1/dt+0.5*mu*(KX2+KY2))*((1/dt+f(wf)+0.5*mu*(KX2+KY2))*wf)
        
        #---double steps with integrating factor method(4th-order Runge-Kutta)---#
        wf=np.exp(-mu*(KX2+KY2)*dt)*wf 
        g1=f(wf          )
        g2=f(wf+0.5*dt*g1)
        g3=f(wf+0.5*dt*g2)
        g4=f(wf+    dt*g3)
        wf=wf+dt*(g1+2*g2+2*g3+g4)/6
                
        if(it%isav==0):
            w   =np.real(sf.ifft2(wf))
            whst[it//isav,:,:]=w         
            
    return locals()

def f(wf):
    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))
    wxf = 1j*KX*wf  ; wx=np.real(sf.ifft2(wxf*KXD*KYD))
    wyf = 1j*KY*wf  ; wy=np.real(sf.ifft2(wyf*KXD*KYD))
    advw =vx*wx+vy*wy
    advwf=sf.fft2(advw)

    
    return -advwf#-mu*(KX2+KY2)*wf
In [3]:
nx=512; ny=512; nt=1200; isav=10
mu=0.0002; dt=0.1
w=np.random.randn(nx,ny)*5

dat=VTE(nx,ny,nt,dt,mu,w,isav)
locals().update(dat)

  0%|          | 0/1200 [00:00<?, ?it/s]/home/mnakanot/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:44: RuntimeWarning: divide by zero encountered in true_divide
/home/mnakanot/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:44: RuntimeWarning: invalid value encountered in true_divide
100%|██████████| 1200/1200 [11:00<00:00,  1.82it/s]
In [4]:
plt.rcParams['animation.html'] = 'html5'
fig=plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

def update_anim(it):    
    ax.clear()

    ax.imshow(whst[it,:,:],origin='lower',aspect='auto',cmap='jet')
    ax.set_title(r"$\omega$")
    ax.set_xlabel(r'$x$')
    ax.set_ylabel(r'$y$')
    
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)    
plt.close()
anim
Out[4]: