【Python4CFD】笔记step1-4

CFDpython - 12 steps to N-S equation

最近在等师兄数据，偶然发现github上有一个很有意思的项目，用python（其实是jupyter）学习CFD，借此机会来摸摸鱼，顺便记录一下做点笔记

00-前言

1. 课程链接CFD Python, a.k.a. the 12 steps to Navier-Stokes
2. 视频链接ME702-Youtube
3. 需要的库环境：

ipywidgets==7.2.1
jupyter==1.0.0
numpy==1.14.3
matplotlib==2.2.2
scipy==1.1.0
sympy==1.1.1

3. 前言部分主要介绍了python+numpy的简单应用，包括以下六个部分：

00-1 Libraries

1. 这一部分就是介绍怎么安装库，个人建议用pycharm+anaconda
2. 注意不要用from xxx import *
3. linspace是一个很有用的函数，用来建立等距的行向量

# np.linspace
myarray = numpy.linspace(0, 5, 10)
# 0-5之间产生10个点
# array([ 0.        ,  0.55555556,  1.11111111,  1.66666667,  2.22222222,
#        2.77777778,  3.33333333,  3.88888889,  4.44444444,  5.        ])

00-2 Variables

1. python无需像c一样声明变量的数据类型
2. 但是要注意，整数/整数=浮点数（int/int=float)

a = 5        #a is an integer 5
b = 'five'   #b is a string of the word 'five'
c = 5.0      #c is a floating point 5  
#%%
type(a) # int
#%%
type(b) # str
#%%
type(c) # float

00-3 Whitespace in Python

1. python中的语句关系由缩进来控制的，而不是像c一样的中括号
2. 这个的话用的例子是loop

00-4 Slicing Arrays

1. 切片，和c一样，从0开始，N-1结束
2. 冒号表示从start到end-1

myvals = numpy.array([1, 2, 3, 4, 5])

myvals[0], myvals[4] # output:1,5

myvals[5] # output: error

myvals[0:3] # output:1 2 3

00-5 Assigning Array Variables

1. 感觉像是python里面比较坑的机制了：引用机制，如果中途修改变量，会导致原变量也进行修改
2. 这个时候需采用copy来解决

a1 = numpy.linspace(1,5,5) # a1=1 2 3 4 5

b = a1

b[3] = 17 # b = a1 = 1 2 3 17 5

a2 = numpy.linspace(1,5,5) # a1=1 2 3 4 5

c = a2.copy()

c[3] = 17 # a2=1 2 3 4 5; c=1 2 3 17 5

00-6 learn more

这一节学到个新东西，原来jupyter里面还可以展示youtube，下次试试能不能展示bilibili

from IPython.display import YouTubeVideo
# a short video about using NumPy arrays, from Enthought
YouTubeVideo('...')

01-Step 1: 1-D Linear Convection 一维线性对流方程

1. 对于一维线性对流方程：
   $$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0$$
2. 其本质是一个双曲型偏微分方程，假设初值$u=u_0$根据特征线可以知道其解析解为:
   $$u=u_0(x\pm ct)$$
3. 该例子是对时间域采用向后一阶差分，对空间域采用向前一阶差分，通过有限差分法对方程进行离散：
   $$\frac{u_i^{n+1}-u_i^n}{\Delta t}+c\frac{u_i^n-u_{i-1}^n}{\Delta x}=0$$
   $$u_i^{n+1}=u_i^n-c\frac{\Delta t}{\Delta x}\left(u_i^n-u_{i-1}^n\right)$$
4. 所有代码如下：

# Remember: comments in python are denoted by the pound sign
import numpy                       #here we load numpy
from matplotlib import pyplot      #here we load matplotlib
import time, sys                   #and load some utilities

#this makes matplotlib plots appear in the notebook (instead of a separate window)
%matplotlib inline   

# 绘制网格：空间节点有41个，时间节点为25个，可以理解为解域是一个41*25的均匀网格
nx = 41  # 加大这个值，会让结果更为精确（即加密网格，但不是所有方程都会精确）
dx = 2 / (nx-1)
nt = 25    #nt is the number of timesteps we want to calculate
dt = .025  #dt is the amount of time each timestep covers (delta t)
c = 1      #assume wavespeed of c = 1

# 设置初始条件：在x属于0.5-1的时候是2，其余时刻是1
u = numpy.ones(nx)      #numpy function ones()
u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s
pyplot.plot(numpy.linspace(0, 2, nx), u);

# 设置temp变量，用来表示上一时刻的值，即n时刻，最终计算n+1时刻的值
un = numpy.ones(nx) #initialize a temporary array

# 开始求解方程:原作者说，这种方式比较浪费资源，之后课程有更好的方式
for n in range(nt):  #loop for values of n from 0 to nt, so it will run nt times
    un = u.copy() ##copy the existing values of u into un
    for i in range(1, nx): 
        u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])

pyplot.plot(numpy.linspace(0, 2, nx), u);

5. 网格可以通过如下代码进行观察：

xx=numpy.linspace(0,2,nx)
tt=numpy.linspace(0,dt*nt,nt)

x,y=numpy.meshgrid(xx,tt,indexing='xy')
pyplot.scatter(x,y)

02-Step 2: Nonlinear Convection非线性对流方程

1. 用u替换原来的c：
   $$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=0$$
2. 其有限差分形式为：
   $$u_i^{n+1}=u_i^n-u_i^n\Delta x$$