【Python4CFD】笔记step5-8

CFDpython - 12 steps to N-S equation

承接上文，这是step5-8，相当于课程的第二阶段，这一阶段主要是讲2维的，先老实把链接放在这：

1. 课程链接CFD Python, a.k.a. the 12 steps to Navier-Stokes
2. 视频链接ME702-Youtube

06 Array Operations with NumPy

1. 在开始前，作者介绍了numpy的数组操作，举了一个很简单的例子，切片计算：

import numpy

u = numpy.array((0, 1, 2, 3, 4, 5))

# 方法一：一个一个算相邻元素的差值
for i in range(1, len(u)):
    print(u[i] - u[i-1])

# 方法二：利用切片直接计算
u[1:] - u[0:-1]
# (array([1, 2, 3, 4, 5])-array([0, 1, 2, 3, 4]))

2. 其次，作者还重点强调了，使用numpy内置的数组和函数是可以加速运算的（原因在于numpy的底层是C内核，并且在数据结构和算法上做了优化）
3. 这里需要介绍一下iPython的一个magic function–>%%timeit
   • 这个函数会返回jupyter一段函数的平均执行时间
   • 注意嗷，这个函数是先运行，再计算，如果在这段section里画图了，会重复画一次图！
4. 例子的代码如下：

#%% 网格划分：求解的是二维对流方程
nx = 81
ny = 81
nt = 100
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = .2
dt = sigma * dx

x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 2, ny)

u = numpy.ones((ny, nx)) ##create a 1xn vector of 1's
un = numpy.ones((ny, nx)) 

### 初始条件

u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

#%% 例子1:raw python
%%timeit
u = numpy.ones((ny, nx))
u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

for n in range(nt + 1): ##loop across number of time steps
    un = u.copy()
    row, col = u.shape
    for j in range(1, row):
        for i in range(1, col):
            u[j, i] = (un[j, i] - (c * dt / dx * 
                                  (un[j, i] - un[j, i - 1])) - 
                                  (c * dt / dy * 
                                   (un[j, i] - un[j - 1, i])))
            u[0, :] = 1
            u[-1, :] = 1
            u[:, 0] = 1
            u[:, -1] = 1


#%% 例子2:python with numpy
%%timeit
u = numpy.ones((ny, nx))
u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

for n in range(nt + 1): ##loop across number of time steps
    un = u.copy()
    u[1:, 1:] = (un[1:, 1:] - (c * dt / dx * (un[1:, 1:] - un[1:, 0:-1])) -
                              (c * dt / dy * (un[1:, 1:] - un[0:-1, 1:])))
    u[0, :] = 1
    u[-1, :] = 1
    u[:, 0] = 1
    u[:, -1] = 1

5. 以上两个例子的结果如下(我的设备是Macbook air M2，不同机子的运行时间不一样的)
   • 例子1:1.63 s ± 31.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
   • 例子2:3.32 ms ± 39.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

07 Step 5: 2-D Linear Convection 2维线性扩散

071 简单的有限差分

1. 这一步很简单，就是把原来的1D变成了现在的2D，方程如下：
   $$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}+c\frac{\partial u}{\partial y}=0$$
   $$\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t}+c\frac{u_{i,j}^n-u_{i-1,j}^n}{\Delta x}+c\frac{u_{i,j}^n-u_{i,j-1}^n}{\Delta y}=0$$
2. 其中，i表示x方向的index，j表示y方向的index，n表示时间上的index
3. 例子给的初始条件是：$u(x,y)=\{\begin{array}{l}\begin{array}{cc}2\text{for}& 0.5\le x,y\le 1\\ 1\text{for}& \text{everywhere else}\end{array}\end{array}$
4. 边界条件：$u=1\text{for }\{\begin{array}{l}\begin{array}{c}x=0,2\\ y=0,2\end{array}\end{array}$
5. 原文给了初始条件的3D图，还是很赞的3D图
6. 代码如下：

from mpl_toolkits.mplot3d import Axes3D    ##New Library required for projected 3d plots

import numpy
from matplotlib import pyplot, cm
%matplotlib inline

###variable declarations
nx = 81
ny = 81
nt = 100
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = .2
dt = sigma * dx

x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 2, ny)

u = numpy.ones((ny, nx)) ##create a 1xn vector of 1's
un = numpy.ones((ny, nx)) ##

###Assign initial conditions
##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2 

###Plot Initial Condition
##the figsize parameter can be used to produce different sized images
fig = pyplot.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')                      
X, Y = numpy.meshgrid(x, y)                            
surf = ax.plot_surface(X, Y, u[:], cmap=cm.viridis)

# array operation
u = numpy.ones((ny, nx))
u[int(.5 / dy):int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

for n in range(nt + 1): ##loop across number of time steps
    un = u.copy()
    u[1:, 1:] = (un[1:, 1:] - (c * dt / dx * (un[1:, 1:] - un[1:, :-1])) -(c * dt / dy * (un[1:, 1:] - un[:-1, 1:])))
    u[0, :] = 1
    u[-1, :] = 1
    u[:, 0] = 1
    u[:, -1]