[問題討論]使用Python學習CFD初級理論系列一數組操作（7/10）

Oler

2019年4月28日 15:13

瀏覽：2721 評論：1

前面的案例中大量采用了Python數值計算包numpy，然而并未使用到numpy的性能。numpy的數值計算實際上調用的是c語言操作，按道理計算速度應該不會慢才對。

numpy的數組操作

在計算量集中的程序中，使用numpy內置的函數操作能夠有效地提高計算性能。下面來舉一個例子，考慮到CFD中經常會遇到如下的迭代式：

[問題討論]使用Python學習CFD初級理論系列一數組操作（7/10）的圖1

假設給定初始值 [問題討論]使用Python學習CFD初級理論系列一數組操作（7/10）的圖2 ，可以通過迭代計算得到的值。

采用迭代方法的代碼可寫成以下形式。

import numpy as np
u = np.array([0,1,2,3,4,5])
un= u.copy()
for i in range(1,len(u)):
   print(u[i] - u[i-1])

輸出結果為：

其實可以改用numpy內置數組操作來實現，代碼寫成以下形式。

import numpy as np
u = np.array([0,1,2,3,4,5])
print(u[1:] - u[0:-1])

輸出結果為：

[1 1 1 1 1]

兩者結果一致。這里采用numpy數組分片功能來進行計算，來看看u[1:]與u[0:-1]到底是多少。

import numpy as np
u = np.array([0,1,2,3,4,5])
print(u[1:])
print(u[0:-1])

輸出結果為：

[1 2 3 4 5]
[0 1 2 3 4]

可以看到u[1:]取的是第2個元素到最后一個元素的值，而u[0:-1]取得的是第1個元素到倒數第2個元素的值。

注：關于數組分片，可參閱numpy官方文檔中的說明。

這樣做有什么用呢？對于大量密集計算來講，這樣做能節省大量的計算時間，下面來用一個復雜案例進行測試。

測試案例

測試案例采用二維穩態導熱問題，其控制方程為拉普拉斯方程：

[問題討論]使用Python學習CFD初級理論系列一數組操作（7/10）的圖4

采用二階中心差分，離散方程可寫成：

[問題討論]使用Python學習CFD初級理論系列一數組操作（7/10）的圖5

很容易得到：

[問題討論]使用Python學習CFD初級理論系列一數組操作（7/10）的圖6

不說二話，直接上代碼。

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

dx = 0.1
dy = 0.1
dx2 = dx*dx
dy2 = dy*dy

def laplace(u):
   nx,ny = u.shape
   for i in range(1,nx-1):
       for j in range(1,ny-1):
           u[i,j] = ((u[i+1,j]+ u[i-1,j]) * dy2 + + (u[i,j+1]+u[i,j-1])* dx2) /(2*(dx2+dy2))


def mat_laplace(u):
   u[1:-1,1:-1]= ((u[2:,1:-1]+u[:-2,1:-1])*dy2 + (u[1:-1,2:] + u[1:-1,:-2])*dx2) / (2*(dx2+dy2))


def calc(N,Niter=100):
   u = np.zeros([N,N]) #全局初始化
   u[0] = 1   #邊界條件
   for i in range(Niter):
       laplace(u)
   return u

def mat_calc(N,Niter = 100):
   u = np.zeros([N,N])
   u[0] = 1
   for i in range(Niter):
       mat_laplace(u)
   return u


# 運行測試
start  = time.clock()
U = calc(100,3000)
end = time.clock()
time_elapsed = end - start #普通計算耗時


start  = time.clock()
U = mat_calc(100,3000)
end = time.clock()
time_elapsed1 = end - start #數組計算耗時


plt.pcolormesh(U)
plt.show()
plt.legend()
print(time_elapsed)
print(time_elapsed1)

計算結果如圖所示。

[問題討論]使用Python學習CFD初級理論系列一數組操作（7/10）的圖7

看到差距了沒有呢，采用循環計算耗時64.14686s，而利用數組運算用時0.44228s，速度提升了一百多倍。這里采用的網格是100*100，迭代次數3000次，如果是真實的CFD計算，網格動不動幾百萬幾千萬，那差距可就海了去了。所以，在利用Python進行數值密集型計算時，盡可能使用numpy內置函數操作，減少循環操作是非常有必要的。關于numpy的具體應用，可上其官網查看。

注：本系列教程來自國外一個使用Python進行CFD初級理論學習的項目，源項目網址為：http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/。感興趣的同學可以去官方主頁了解更多信息。

本文轉載自微信公眾號“CFD之道”，有刪減，感謝源作者。

登錄后免費查看全文

立即登錄