Visualize cavity flow with matplotlib and save as gif animation

http://nbviewer.ipython.org/github/barbagroup/CFDPython/blob/master/lessons/15_Step_11.ipynb

I found out how to write the cavity flow, but it's not enough because it doesn't animate.

Below is the code that displays only the final result. (In the original story, I use iPython Notebook)

cavorg.py

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation

nx = 41
ny = 41
nt = 5
nit=50
c = 1
dx = 2.0/(nx-1)
dy = 2.0/(ny-1)
x = np.linspace(0,2,nx)
y = np.linspace(0,2,ny)
X,Y = np.meshgrid(x,y)

rho = 1
nu = .1
dt = .001

u = np.zeros((ny, nx))
v = np.zeros((ny, nx))
p = np.zeros((ny, nx)) 
b = np.zeros((ny, nx))


def buildUpB(b, rho, dt, u, v, dx, dy):
    
    b[1:-1,1:-1]=rho*(1/dt*((u[1:-1,2:]-u[1:-1,0:-2])/(2*dx)+(v[2:,1:-1]-v[0:-2,1:-1])/(2*dy))-\
                      ((u[1:-1,2:]-u[1:-1,0:-2])/(2*dx))**2-\
                      2*((u[2:,1:-1]-u[0:-2,1:-1])/(2*dy)*(v[1:-1,2:]-v[1:-1,0:-2])/(2*dx))-\
                      ((v[2:,1:-1]-v[0:-2,1:-1])/(2*dy))**2)

    return b
    
def presPoisson(p, dx, dy, b):
    pn = np.empty_like(p)
    pn = p.copy()
    
    for q in range(nit):
        pn = p.copy()
        p[1:-1,1:-1] = ((pn[1:-1,2:]+pn[1:-1,0:-2])*dy**2+(pn[2:,1:-1]+pn[0:-2,1:-1])*dx**2)/\
                        (2*(dx**2+dy**2)) -\
                        dx**2*dy**2/(2*(dx**2+dy**2))*b[1:-1,1:-1]

        p[-1,:] =p[-2,:] ##dp/dy = 0 at y = 2
        p[0,:] = p[1,:]  ##dp/dy = 0 at y = 0
        p[:,0]=p[:,1]    ##dp/dx = 0 at x = 0
        p[:,-1]=0        ##p = 0 at x = 2
        
    return p
    
def cavityFlow(nt, u, v, dt, dx, dy, p, rho, nu):
    un = np.empty_like(u)
    vn = np.empty_like(v)
    b = np.zeros((ny, nx))
    
    for n in range(nt):
        un = u.copy()
        vn = v.copy()
        
        b = buildUpB(b, rho, dt, u, v, dx, dy)
        p = presPoisson(p, dx, dy, b)
        
        u[1:-1,1:-1] = un[1:-1,1:-1]-\
                        un[1:-1,1:-1]*dt/dx*(un[1:-1,1:-1]-un[1:-1,0:-2])-\
                        vn[1:-1,1:-1]*dt/dy*(un[1:-1,1:-1]-un[0:-2,1:-1])-\
                        dt/(2*rho*dx)*(p[1:-1,2:]-p[1:-1,0:-2])+\
                        nu*(dt/dx**2*(un[1:-1,2:]-2*un[1:-1,1:-1]+un[1:-1,0:-2])+\
                        dt/dy**2*(un[2:,1:-1]-2*un[1:-1,1:-1]+un[0:-2,1:-1]))

        v[1:-1,1:-1] = vn[1:-1,1:-1]-\
                        un[1:-1,1:-1]*dt/dx*(vn[1:-1,1:-1]-vn[1:-1,0:-2])-\
                        vn[1:-1,1:-1]*dt/dy*(vn[1:-1,1:-1]-vn[0:-2,1:-1])-\
                        dt/(2*rho*dy)*(p[2:,1:-1]-p[0:-2,1:-1])+\
                        nu*(dt/dx**2*(vn[1:-1,2:]-2*vn[1:-1,1:-1]+vn[1:-1,0:-2])+\
                        (dt/dy**2*(vn[2:,1:-1]-2*vn[1:-1,1:-1]+vn[0:-2,1:-1])))

        u[0,:] = 0
        u[:,0] = 0
        u[:,-1] = 0
        u[-1,:] = 1    #set velocity on cavity lid equal to 1
        v[0,:] = 0
        v[-1,:]=0
        v[:,0] = 0
        v[:,-1] = 0
   
    return u, v, p
    
fig = plt.figure(figsize=(11,7), dpi=100)
u, v, p = cavityFlow(nt, u, v, dt, dx, dy, p, rho, nu)

plt.contourf(X,Y,p,alpha=0.5)
plt.colorbar()
plt.contour(X,Y,p)
plt.quiver(X[::2,::2],Y[::2,::2],u[::2,::2],v[::2,::2])
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

It seems that you can make a gif animation using matplotlib.animation, so I tried it, but it didn't work. The reason will be explained below.

You can make gif animation if imagemagick is included. See below for this. http://kiito.hatenablog.com/entry/2013/11/27/233507

The specific method for animating is as follows.

There seem to be two ways to animate. animation.ArtistAnimation ... Draw the data prepared in advance animation.FuncAnimation ... Update data from time to time

from http://cflat-inc.hatenablog.com/entry/2014/03/17/214719

Other examples http://qiita.com/HirofumiYashima/items/5864fb9384053b77cf6f http://d.hatena.ne.jp/xef/20120629/p2 http://stackoverflow.com/questions/24166236/add-text-to-image-animated-with-matplotlib-artistanimation

Since the first simulation code is running a for loop in cavityFlow (), I tried to spit out data for each loop and draw with ArtistAnimation, but when I pass contourf to the array of drawing data, it is a bug. coming out. This has been reported elsewhere, but I'm not sure if it's supported.

http://matplotlib.1069221.n5.nabble.com/Matplotlib-1-1-0-animation-vs-contour-plots-td18703.html http://stackoverflow.com/questions/23070305/how-can-i-make-an-animation-with-contourf

For the time being, if you do it with Func Animation instead of Artist Animation, it seems that contourf will also be an animation, but there is also information that you can do it with Artist Animation. Except for contourf, I tried to draw only streamlines and it was created successfully, so this seems to be a bug caused by contourf.