面向邊緣計算的神經網絡MPS加速算法*

王卓霖， 江俊揚， 楊耿超， 姚清河， 蔣子超， 張儀

1.中山大學航空航天學院，廣東 深圳 518107

2.特文特大學機器人和機電一體化系，荷蘭 恩斯赫德 7522NB

計算流體力學（CFD，computational fluid dynamics）是當前流體力學領域中十分重要的研究方法，通過數值求解流體運動的控制方程組，研究流體流動的物理現象.目前主流計算流體方法主要分為兩類：基于網格的歐拉方法和基于粒子的拉格朗日方法.歐拉方法如有限體積法（FVM）、有限元法（FEM）、有限差分法（FDM）通過對計算區域劃分網格，離散微分方程，獲得網格點各時刻的物理量變化；拉格朗日方法則研究粒子的物理量隨時間的變化，適用于比較復雜的邊界情況（Chikazawa et al.，2001）或者大變形流體（Gotoh et al.， 2006）模擬場景.

近年來拉格朗日方法得到發展，SPH（Monaghan，1992）、MPS（Koshizuka et al.，1996）等粒子方法逐漸受到CFD 界重視.其中MPS（moving particle semi-implicit ）最早由Koshizuka et al.于1996 年提出，并由于其準確性和穩定性等特點受到了學術界的廣泛關注.隨著MPS理論的進一步發展，Khayyer et al.（2008）提出了CMPS方法修正了梯度算子導致的動量不守恒，增加了自由表面粒子的穩定性，通過改進PPE方程和修正矩陣提出了MPS-HS 方法、MPS-HL 方法、ECS 方法、GC 方法等，降低了計算的不穩定性；Tsuruta et al.（2013）提出了MPS-DS 方法修正了斥力模型，進一步增強了計算的穩定性；Tsuruta et al.（2015）提出了SPP 方法通過引入空間勢能粒子對非物理空腔問題做了優化；越來越多的修正方法強化了MPS 的優越性，抑制了壓力振蕩，使其被大量應用于涉及復雜邊界、大變形、多相流等復雜流體問題中（Shibata et al.，2007；張雨新等，2012；楊超，2018）.

MPS 采用預估和修正（顯式和隱式）兩個步驟來求解流體控制方程（Koshizuka et al.，1996），預估步中通過體積力與黏性力顯式更新速度和位移；在校正步中通過保持粒子數密度不變構造壓力泊松方程，隱式求解下一時間步的壓力，并利用壓力梯度修正速度和位移.然而，這種隱式計算方法需消耗大量計算資源.與顯式計算壓力狀態方程的SPH 方法相比，傳統MPS 方法難以高效并行求解，使其不適合大規模仿真計算，在實際工程中無法廣泛應用（張弛等，2011）.前期研究中，我們提出了一種基于數據驅動求解MPS 壓力泊松方程的NN-MPS 算法，該算法將原始MPS 中壓力泊松方程的求解重新描述為回歸問題，提取對泊松求解起主導作用的特征參數作為回歸問題的輸入，并引入神經網絡（NN，neural network）解決該回歸問題.

邊緣計算（何騰，2020）是相對于云計算（崔勇等，2017）而言的一種概念，指在靠近數據源頭的網絡邊緣側對數據進行直接計算，而無需將數據傳輸至云端數據處理中心的一種數據處理方式.隨著人工智能的發展，在數據源頭的離線設備上直接進行實時神經網絡推理的需求增加，昂貴、大體積、大功率的傳統計算硬件如專業級GPU（Lee et al.，2010）、CPU 集群已經無法滿足目前的需求，廉價、高性能的面向神經網絡的邊緣計算硬件成為研究熱點.因此，多種專供神經網絡的專用集成電路（ASIC）芯片被研發出來，如寒武紀思元系列（Liu et al.，2015）、海思科技的昇騰910、谷歌的TPU（Jouppi et al.，2017）及本文將使用的Atlas 200 DK 搭載的Ascend 310 等都屬于這類ASIC 芯片.在工業生產領域，常有在邊緣側對監測數據進行離線分析的場景（如地質災害監測（杜年春等，2022）、空氣質量監測（滕云豪，2022）），而由于邊緣端算力限制，目前在邊緣側利用流體仿真進行數據分析的應用仍較為少見.我們在NN-MPS中對神經網絡的引入，不僅提高了在傳統硬件上的計算速度，也為MPS 方法在這些前沿的邊緣計算硬件上的計算提供了可能.

本文使用電功耗、成本、體積等各方面較有優勢，且在張量計算等神經網絡計算有出色計算效率的邊緣計算設備Atlas 200 DK 作為加速硬件，成功實現了以新型邊緣計算硬件作為加速設備的計算流體仿真.本文研究用較低的成本在小規模與大規模粒子法仿真中達到了與傳統CPU計算效率相似甚至有所提高的計算效果.該研究為將來特殊場景下的快速流體仿真提供了一種可行的參考，意味著我們可以在一些需要低成本或無法使用專業計算集群或是在實時采集實際物理邊界信息的邊緣側場景中，提供一定的快速流體動力學分析能力.

1 移動粒子半隱式（MPS）方法

1.1 控制方程

不可壓縮流體的控制方程組由連續性方程和Navier-Stokes方程組成：

其中ρ為流體密度，P為壓力，v為速度向量，f為流體單位質量力，ν是運動黏度系數.

1.2 粒子作用模型

1.2.1 核函數MPS方法中，核函數與光滑半徑內鄰域粒子相互作用，計算控制方程中的微分算子.Koshizuka et al.（1996）提出的核函數由如下表示：

式中re為粒子的有效半徑，即光滑長度.

核函數值在靠近中心處最大，并隨著r的增加減小.因此，核函數在計算中起到了類似于權重函數的作用，使較遠粒子對于當前粒子的影響較小.在本節后面的部分中將展示由權重函數表示的微分離散計算方法.

1.2.2 粒子密度模型MPS 中粒子數密度n i是第i粒子在核函數控制域內與其周圍粒子核函數數值的累加：

其中N為粒子i的相鄰粒子，W(|ri-rj|，h)為核函數，ri、rj分別為i、j粒子的位置矢量，h為光滑長度.粒子數密度與流體密度成正比.

1.2.3 Laplacian模型Laplacian算子采用如下形式：

式中D為計算維數，

1.3 預估-修正法

為求解由式（1）～（2）組成的方程組，MPS 方法通過引入中間速度，將N-S 方程分兩步求解.第一步通過外力信息和流體速度信息，顯式求解中間速度v*，第二步通過壓力梯度修正速度，即

其中

1.4 壓力泊松方程（Pressure Poisson Equation）

傳統MPS 算法通過保持粒子數密度不變來實現流體的不可壓縮性，利用中間密度ρ*和中間速度v*，令下一時間步的速度散度為0、密度變化為0，可構造壓力泊松方程，進而求解Pk+1.參考Tanaka et al.（2010）、Lee et al.（2011）等的做法，本文采用基于散度的表達式：

其中右邊第一項為基于散度為0 構造泊松方程的源項、第二項為基于密度不變構造的泊松方程源項.γ為可變參數，本文中定義為0.2.

1.5 時間積分

綜上，MPS算法每個時間步的求解步驟如下：

（i）根據式（8）中黏性力和質量力獲得中間速度v*，并計算中間位置r*；

（ii）根據式（4）計算中間粒子數密度n*；

（iii）根據方程（9）求解壓力泊松方程，得到Pk+1；

（iv）利用所求壓力梯度獲得修正速度v'，并計算vk+1、rk+1.

2 基于神經網絡的MPS算法（NN-MPS）

2.1 原始壓力泊松方程

根據式（5）離散Laplacian模型，PPE左邊項可以離散成如下形式：

代入式（9），并改成矩陣形式可以表示為

其中

其中下標n為參與計算的粒子數.可知，PPE 是一個包含所有粒子狀態的矩陣方程，在大規模計算的情況下，求解該方程將需要非常大的存儲空間和非常長的求解時間.為了更高效地求解PPE、MPS中常用共軛梯度法（CG）、不完全Cholesky共軛梯度法（ICCG）進行求解大型矩陣方程，求解該方程的代價仍然極高.

2.2 數據驅動的壓力泊松方程

為了解決PPE 求解代價過高的問題，前期研究中，我們提出了一種數值驅動的NN-MPS 方法.NNMPS方法將求解PPE方程轉化為機器學習領域的回歸問題.針對這一回歸問題，算法采用神經網絡建立求解模型.對于高維問題，傳統數值方法的計算成本過高，神經網絡較強的提取特征的能力以及求解非線性問題的能力，使其成為了解決該回歸問題一個很好的選擇.

為了使回歸求解器更加準確和穩定，我們需要為模型選擇合適的輸入.參考原始MPS 算法的PPE構造：

觀察式（15），可以發現中心粒子的壓力依賴于光滑半徑內其他粒子的壓力權重、自身的中間速度、中間粒子數密度.據此，NN-MPS算法提出一種回歸模型的構造形式：

其中

式中Pk表示當前時間步該粒子的壓力，則體現了周圍粒子與中心粒子的壓力積分特征，的數量由模型的鄰域極限半徑確定.

基于以上數據模型，NN-MPS構造了適合PPE回歸求解器的輸入.

2.3 DNN（Deep Neural Network）網絡模型

多層感知機（Multi-Layer Perceptron）作為一種經典深度學習模型，是當前主流深度神經網絡的基礎（Gardner et al.，1998）.一個多層感知機包含一個輸入層、一個輸出層，以及多個隱藏層，每一層包括多個感知器，每個感知器具有一個或多個輸入、偏置、激活函數，以及單個輸出.感知機每一層可以表示成：

其中x=(x1，x2，…，xn)表示該層的輸入向量，ω表示權重，b表示偏置，H則稱為激活函數，f(x)表示該層輸出.

本文選擇以多層感知機作為回歸網絡的基礎結構.隱藏層的層數對回歸求解效果、計算速度有較大的影響，目前相關領域尚無法給出DNN 模型深度（隱藏層層數）與計算誤差之間的絕對關系.因此，綜合考慮了計算精度、計算效率之間的平衡，結合實際數值測試結果，我們選擇將隱藏層數設置為5.如圖1 所示，整個結構包含1 個輸入層、5 個隱藏層和1 個輸出層，層與層之間采用ReLU（Glorot et al.，2011）函數作為激活函數.輸入層中的排序由粒子與中心粒子距離決定，這種排序方式可以幫助神經網絡學習系統的潛在特征.

圖1 神經網絡結構示意圖Fig.1 Structure of the neural network

改進后的MPS算法單時間步積分流程如下：

（i）根據式（8）中黏性力和質量力獲得中間速度v*，并計算中間位置r*；

（ii）根據式（4）計算中間粒子數密度n*；

（iii）構造神經網絡輸入，并通過神經網絡求解Pk+1；

（iv）利用求得的壓力，根據式（9）獲得修正速度v'，并計算vk+1，rk+1.

3 利用邊緣計算設備加速MPS算法

3.1 Atlas 200 DK開發板

華為公司的Atlas 200 DK 開發者套件是以Atlas 200 AI 加速模塊為核心的開發板形態.其中Atlas 200 AI 加速模塊集成了Ascend 310 AI 處理器，Ascend 310 芯片內置2 個DaVinci AI core，可以實現最大8TFLOPS/FP16或16TOPS/INT8的乘加性能（華為技術有限公司，2022）.

開發板可以實現圖像、視頻等多種數據分析與推理計算，目前廣泛用于基于機器視覺的AI應用，如智能監控、機器人、無人機、視頻服務器等場景.機身包括2 個DaVinci AI core，1 個8 核A55CPU，支持1000 M 以太網.由于其對神經網絡的高速計算的特性與高速以太網通信能力，我們選擇Atlas 200 DK作為配置NN-MPS的邊緣計算設備.

本文所有工作基于Atlas 200 DK開發者套件的1.32.0.0版本進行開發.

3.2 開發流程

本文采用Python3在Mind Studio平臺上進行NN-MPS程序開發.開發前開發平臺為Mind Studio，語言為Python3、hiai庫。

華為公司為調用AI Core 提供了一套如圖2 所示的昇騰軟件棧，其相關概念的具體介紹可見華為官方文檔.在利用其中的模型轉換工具將Tensorflow 模型轉換為其所需格式后，通過如下的流程實現調用AI Core完成神經網絡推理：

圖2 昇騰AI軟件棧Fig.2 Software stack of Ascend AI

① 初始化圖（Graph）對象，圖的作用是管理和編排計算引擎；

② 創建引擎（Engine），引擎是執行功能的基本單元；

③ 設置模型文件（.om）路徑；

④ 構造輸出輸入對象；

⑤ 調用引擎進行模型推理；

⑥ 結果后處理.

3.3 Socket通信方案

Atlas開發板板載CPU實際性能并不比一般計算機CPU高，因此將整個NN-MPS由開發板運行并不是個很合適的方案.考慮到Atlas 200 DK支持1000 M以太網傳輸，將計算過程分成兩個部分，推理部分由開發板完成，其余部分由電腦CPU 完成具有相當高的可行性.因此，為了結合計算機CPU 算力與邊緣計算設備對神經網絡的加速能力，程序將求解壓力泊松方程的神經網絡推理部分在Atlas 200 DK 開發板上實現，其余部分由計算機CPU 完成，兩者通過千兆以太網利用Socket通信協議實現快速的數據交換.圖3展示了整體程序結構.

圖3 整體計算流程Fig.3 Calculation flow chart

4 算例實驗

4.1 模型描述

本文中采用如圖4所示增加障礙物的潰壩模型作為驗證算例，模型由傳統潰壩模型加上圓柱形障礙物組成.文中共有3種模型：潰壩加固定圓柱障礙物、潰壩加周期性移動圓柱障礙物，以及增大寬度的潰壩加固定圓柱障礙物.其中，移動圓柱障礙物的模型則是在圖4（a）的基礎上為圓柱增加vx= 2.5*2*sign((step + 10)%19 - 10)的運動率，式中step 為計算時間步，sign(a)代表取a的符號（+/-）.神經網絡用的訓練集為傳統MPS 計算經典無障礙物的潰壩模型過程中輸出每時間步的壓力積分特征、速度散度、粒子數密度與中心粒子下一時間步壓力所得.

圖4 MPS模型初始狀態Fig.4 MPS model initial state

采用該模型的原因主要有以下幾點：首先，潰壩模型為大部分MPS 算例使用的經典模型，且模型結構簡單易于實現；其次，增加了圓柱障礙物和移動障礙物可以更好表現MPS 算法對于復雜邊界條件的適應性和泛化能力；不同的粒子規模有助于比較不同規模流場場景下計算速度的差別.

4.2 誤差分析

本文計算了在小模型場景下，固定障礙物與移動障礙物模型計算結果在經典MPS 算法、使用CPU 的NN-MPS 算法、使用Atlas 200 DK 加速的NN-MPS 算法的計算精度.使用如下3 個參數進行數值比較：近壁面液面高度、液面前緣位置，以及如圖5所示區域的平均壓力.

圖6～8是3個參數在兩種模型下的計算結果，其中前緣位置與液面高度只取了流體與前壁碰撞前的數據.比較兩個NN-MPS 與傳統MPS 的數據，使用NN-MPS 算法在液面高度、前緣位置上有最大在0.1 作有的相對誤差，在壓力結果上雖然與原始MPS 算法相比產生了漂移與變形，相對誤差較大，但是在總體變化趨勢上與實際結果仍比較相似，這些差別是由神經網絡模型的結構或訓練參數導致，可以通過優化模型得到改善；在兩種硬件運行的NN-MPS算法計算結果之間也有微小的差別，這是因為Atlas上的AI芯片只支持最大FP16浮點數，而CPU上則支持64位浮點運算，整體來看，由硬件區別導致的誤差與由神經網絡模型導致的誤差相比，前者幾乎可以忽略.

圖6 前緣位置變化Fig.6 Movement of the leading edge

圖7 近壁面液面高度變化Fig.7 The liquid level height near the left wall

圖8 壓力變化Fig.8 Pressure developments

單獨分析在不同硬件上的NN-MPS計算結果，發現計算前期階段在CPU 與Atlas上的NN-MPS計算結果差別較小，而在第75時間步左右移動圓柱下的前緣位置參數Atlas 與CPU 計算結果開始出現差別；在125 時間步附近，兩者在靜止圓柱和移動圓柱的液面高度也開始出現差別.這可能是由于經過碰撞之后產生的流場變化加劇，導致由FP64移植到FP16設備上產生微小誤差的積累增大，從而使結果出現較大差距.這點可以通過靜止圓柱與移動圓柱計算結果的區別加以證明：靜止圓柱場景下，流場變化與移動圓柱場景相比較為緩和，因而移動圓柱場景下，兩種NN-MPS 算法的計算結果更早出現較大差別，同理在移動圓柱場景下，液面高度也比靜止圓柱場景下出現了更大的差距.

4.3 速度對比

本文比較了加寬模型與小模型分別采用傳統MPS、神經網絡改進后的MPS、移植到Atlas 200 DK的MPS在求解PPE時的用時情況.為了時間比較的合理性，移植到Atlas 的MPS 計算PPE 求解時包括了Socket通信消耗的時間.

實例計算使用的CPU：Intel（R） Core（TM）i7-10510U CPU@1.80 GHz.

表1 展示了神經網絡改進的MPS 算法與移植到Atlas 200 DK 上的改進的MPS 算法在求解PPE的時間上都相較傳統MPS 算法有10 倍以上的提升，數據規模較小的模型中（如算例1、算例2中），Socket 通信耗時占比較大，通信時間與AI 芯片神經網絡求解時間處于同一量級，因而在CPU 與在Atlas上面的MPS-NN算法相比，反而CPU的計算速度較快；而在數據規模較大的模型中（算例3），Socket通信時間占比減小，移植到Atlas 200DK上的MPS 相較在CPU 上的改進MPS 算法求解PPE用時又有了16%的提升.

表1 PPE求解時間比較Table 1 Computational time for solving PPE

5 總 結

本文采用神經網絡結合邊緣計算硬件的方式，將MPS 方法中的計算瓶頸——壓力泊松方程的計算時間縮減，提出了一種面向邊緣計算場景的流體動力學仿真新方案.研究利用神經網絡抽象壓力泊松方程的求解過程，將其轉化為光滑半徑內粒子特征與下一時間步壓力值之間的回歸問題，成功減少了泊松方程求解時間.同時，本文利用邊緣計算硬件求解神經網絡推理過程，探究在邊緣計算場景下進一步加速計算的可能性.在實驗結果中，當粒子規模在28 800 以上時，使用邊緣計算硬件的計算速度相對于使用CPU 的計算速度有所提升，而在小規模場景下，使用邊緣計算硬件與CPU 計算速度相當，因此在計算規模較大時，可以使用本文所述方法作為加速手段.本文采用的方案，相較于使用GPU 等常見加速硬件，具有更低的功率需求和更低的成本，為一些邊緣側有流體動力學仿真需求的場景提供了一種有效的解決方案，同時也為邊緣計算與基礎學科進行學科交融提供新的切入方向.