基于ARM SVE的光滑粒子流體動力學SIMD加速方法*

2021-06-25 09:46:04范小康夏澤宇龍思凡楊燦群

計算機工程與科學 2021年6期

關鍵詞：方法

范小康，夏澤宇，龍思凡，楊燦群

(國防科技大學計算機學院,湖南長沙 410073)

1 引言

數值模擬技術是工程應用和科學研究的一個重要手段。與實驗手段相比，數值模擬成本低，并且可以模擬很多難以開展實驗或無法開展實驗的情況，因而受到越來越多研究人員的青睞。光滑粒子流體動力學SPH(Smoothed Particle Hydrodynamics)[1 - 3]方法是一種無網格的粒子方法，近年來逐漸成為數值模擬技術領域的一個研究熱點。

傳統的數值模擬技術以網格作為基礎，將實際問題轉化為數學描述的離散形式，即物理控制方程組，然后采用計算機求解。比較具有代表性的方法包括有限差分法FDM(Finite Differential Method)[4]、有限體積法FVM(Finite Volume Method)[5]和有限元法FEM(Finite Element Method)[6]，已經在計算流體力學CFD(Computational Fluid Dynamics)和計算固體力學CSM(Computational Solid Mechanics)中取得了廣泛應用。傳統基于網格的方法存在著固有的缺陷，復雜幾何形狀的網格構造是一個難點，此外，基于網格的方法也難以處理大變形問題。

近年來，無網格方法受到了越來越多的重視。無網格方法的主要思想是通過使用一系列任意分布的節點或者粒子來求解具有各種各樣邊界條件的積分方程或偏微分方程組。無網格方法的優點在于免去了網格生成的麻煩，并且避免了網格扭曲與網格重構的問題。SPH是一種無網格的拉格朗日粒子方法，該方法的基本思想是將連續的流體(或固體)用相互作用的質點組來描述，各個質點上承載各種物理量，包括質量、位置坐標、速度等，通過求解質點組的動力學方程和跟蹤每個質點的運動軌道，求得整個系統的力學行為。SPH方法的粒子特性使得其在處理大變形、運動物質表面以及自由表面等問題時優勢明顯，已被諸多商用軟件采用。

作為一種粒子方法，SPH方法在每一個迭代步都需要計算大量的粒子間相互作用，如何提高SPH方法的計算效率是目前的研究難點，也具有非常大的實際意義。針對SPH方法的計算效率問題，國內外也開展了大量的研究工作[7 - 11]。作為一種顯式求解的粒子方法，SPH天然地適合采用并行結構求解，目前已涌現出大量基于分布式存儲的高性能計算平臺、共享存儲結構的并行SPH計算方法，有效地提高了SPH方法的計算效率，使得SPH可處理的問題規模越來越大。然而，以上方法研究的都是SPH方法中的進程級或線程級并行，卻鮮有工作探索SPH方法中的數據級并行。

近年來ARM面向高性能計算推出了可伸縮矢量擴展SVE(Scalable Vector Extension)[12]。SVE并不是ARM NEON的簡單擴展，而是基于ARM AArch64架構的下一代SIMD指令集。與NEON相比，SVE可以獲得超過7倍的加速比，極大地提升了性能[12]。SVE允許硬件設計師根據應用需求定義向量寄存器寬度，最大可支持2 048位的向量寬度。SVE支持VLA(Vector Length Agnostic)編程方式，編程人員或編譯器無需在程序中指定向量寬度，程序被編譯成SVE指令后，可在不同向量寬度的CPU上兼容運行，而無需針對不同向量寬度的CPU重新編譯。SVE在高性能計算領域擁有非常好的應用前景，比如日本的下一代超級計算機富岳[13,14]就采用了SVE。

SVE引入了一些新的架構特點，包括：(1)可變向量長度：SVE提供32個向量寄存器，向量寄存器長度為128的整數倍，最低128位，最高可支持2 048位；(2)每通道預測(Per-Lane Predication):SVE提供16個預測寄存器(Predicate Register)，預測寄存器每一位控制向量寄存器的一個字節，通過預測寄存器中每一位的有效狀態來控制向量寄存器中對應元素是否參與運算；(3)聚集取和分散存(Gather-Load，Scatter Store):支持非連續存儲數據的高效訪問。

SVE提供了比傳統SIMD架構更長的向量寄存器，具備更強的計算能力，在SVE上運行的程序理論上可獲得更高的性能。本文基于SVE提出了一種SPH的SIMD向量優化方法，利用SVE對SPH方法的關鍵部分粒子間相互作用的計算開展向量化并行，探索SPH的數據級并行。實驗結果表明，采用SVE可以獲得粒子間相互作用計算的顯著性能提升。

2 背景介紹

SVE是ARM在AArch64架構基礎上針對高性能應用所推出的SIMD擴展指令集。ARM為SVE定義了32個向量寄存器和16個預測寄存器。向量寄存器的寬度是128的整數倍，范圍從128位～2 048位。預測寄存器的每一位(bit)對應向量寄存器的一個字節(byte)，通過不同位的有效性控制向量寄存器中不同元素是否參與向量操作，預測寄存器可以看作是一個控制向量運算的掩碼。ARM為SVE引入了VLA編程方式，對于循環的向量化可以實現與處理器向量寄存器寬度無關的控制。

2.1 每通道預測

圖1通過一條向量化融合乘加指令(fmla)簡要說明基于預測寄存器的向量操作。指令fmla的操作是把向量寄存器z1與z2中的數據相乘，乘積與向量寄存器z0中的數據相加，最后將結果寫入向量寄存器z0中。向量寄存器的后綴·d表示寄存器中存放的是雙精度浮點(double)類型的操作數。預測寄存器p0/m為合并模式，表示若預測寄存器該有效位無效，則目標寄存器對應通道數據保持原值。3個向量寄存器z0,z1和z2 中參與運算的元素由預測寄存器p0對應的位控制，若對應的位為TRUE(1)，則向量寄存器中對應的元素參與運算，如A[7]、A[5]、A[3]、A[1]，反之若預測寄存器中對應的位為FALSE(0)，則向量寄存器中對應的元素不參與運算，如A[6]、A[4]、A[2]、A[0]。

Figure 1 SIMD operation based on predicate register

2.2 VLA

以DAXPY運算為例簡要說明VLA編程方式。基于VLA的DAXPY運算代碼如下所示：

1 //-------------------------------------------------------------

2 // voiddaxpy(int64_tn,doubleda,double *dx,double *dy) {

3 // for (int64_ti=0;i

4 //dy[i]=dx[i] *da+dy[i];

5 // }

6 // }

7 //-------------------------------------------------------------

8 //x0=n,x1=&dy[0],x2=&dx[0],d0=da

9 daxpy:

10 movx3,0//x3=i=0

11 movz0.d,d0//z0=broadcast(da)

12whileltp0.d,xzr,x0//p0=while(0

13 ·L5:

14 ld1dz1.d,p0/z,[x1,x3,lsl 3]//p0?:z1=dy[i]

15 ld1dz2.d,p0/z,[x2,x3,lsl 3]//p0?:z2=dx[i]

16 fmadz1.d,p0/m,z0.d,z2.d/*p0?:z1 +=dx[i]*da*/

17 st1dz1.d,p0,[x2,x3,lsl 3]//p0?:dy[i]=z1

18 incdx3//i+=(VL/64)

19whileltp0.d,x3,x0//p0=while(i++

20 b.any.L5// more to do?

21ret

代碼中第2～6行顯示的是DAXPY運算的源代碼，第9～21行顯示的是經過編譯之后生成的含有SVE指令的匯編代碼。記處理器的向量寄存器寬度為VL，1次向量操作最多可處理的數據個數為VL/64。第10行和第11行初始化寄存器x3和z0。第12行的“whilelt”指令對預測寄存器p0初始化，分別對比0，1,…,j,…,VL/64-1和n的大小并為p0的對應位賦值，即p0[j]=j

3 SVE向量化優化方法

SPH將計算域內的連續介質離散成一系列具有質量、速度和能量的“粒子”，粒子的物理量通過對其支持域內其他粒子進行插值近似，并將函數導數轉化為一個核函數的導數，從而實現對Navier-Stokes方程組的離散。粒子隨著流體流動并與鄰近的粒子產生交互。SPH方法的計算過程被分解成一系列時間步，每一個時間步更新粒子的物理信息，其中最重要的計算過程是粒子間相互作用的計算。粒子間相互作用計算的串行算法如算法1所示。

算法1SPH粒子間相互作用計算串行算法

Procedure ParticlesInteractionSerial

Begin

1 // 記P為所有粒子的集合

2foreachp∈Pdo

3q←NeighborCellParticles(p);

4for(i=0;i

5d=distance(p,q[i]);

6if(d

7interact(p,q[i]);

8endfor

9endforeach

end

對于計算域P中的任意一個粒子p，首先通過函數NeighborCellParticles獲取該粒子鄰近單元中所有粒子的列表q，對于列表q中的任意一個粒子q[j]，根據p與q[j]的距離d判斷q[j]是否在粒子p的支持域內部。若粒子q[j]在粒子p的支持域內(d

函數interact通過如式(1)～式(3)所示的動量方程、連續性方程和狀態方程計算粒子間相互作用。

動量方程、連續性方程和狀態方程如式(1)～式(3)所示：

(1)

(2)

(3)

3.1 SPH粒子間相互作用向量化并行計算框架

在粒子間相互作用計算的串行算法中，內層循環的每一次迭代都需要判斷2個粒子之間的距離是否小于支持域的半徑，即循環內存在條件分支。條件分支是傳統SIMD優化方法的一個難點，而SVE所提供的預測寄存器則能夠以類似掩碼的形式很好地處理條件分支。

算法2描述了SPH粒子間相互作用計算的SVE向量化算法框架。

算法2SPH粒子間相互作用計算SVE向量化算法

Procedure ParticlesInteractionSVE

Begin

1 // 記P為所有粒子的集合

2foreachp∈Pdo

3q←NeighborCellParticles(p);

4for(i=0;i

5pd0=svwhilelt(i,len(q));

6vq[0:n-1]=ld(pd0,q[i]);

7vd[0:n-1]=distance(pd0,p,vq[0:n-1]);

8pd1=svcmplt(pd0,vd,h);

9interactSVE(pd1,p,vq[0:n-1]);

10endfor

11endforeach

end

如第4行所示，循環的每一次迭代最多處理n個粒子，n的計算方法為向量寄存器的寬度除以單個數據的大小(假如向量寄存器寬度為128位，數據類型為雙精度浮點，即一個數據的大小為64位，則n=128/64=2)。第5行通過循環控制變量j和列表q的長度len(q)構造預測寄存器pd0，預測寄存器pd0包含n個元素，pd0[k]=(j+k)

觀察粒子間相互作用計算的3個方程，動量方程和連續性方程均涉及大量的求和運算，占據了粒子間相互作用的絕大部分計算時間，而且2個方程中存在大量的數據級并行，因此，下面將重點描述如何采用SVE對動量方程和連續性方程進行向量化。以下的向量化并行計算均在預測寄存器的控制下進行。

3.2 動量方程向量化

動量方程是關于粒子速度矢量vk的方程，包含x,y,z3個方向的分量，各對應一個方程，如式(4)～式(6)所示：

(4)

(5)

(6)

Figure 2 SIMD calculation of

Figure 3 SIMD calculation of rij·vij

Figure 4 SIMD calculation of

記:

Figure 5 SIMD calculation of

3.3 連續性方程向量化

Figure 6 SIMD calculation of wij

4 實驗結果

本節對本文所提出的向量化方法的加速效果進行評估。本文采用5個測試用例，分單精度和雙精度2種模式，對比了從128位～2 048位共16種向量寄存器寬度情況下粒子間相互作用向量化算法對串行算法的加速比，得到了單精度和雙精度2種模式下最高3.71倍和4.92倍的加速比。

(1)模擬器。本文采用了ARM指令模擬器ARMIE(ARM Instruction Emulator)對SVE指令開展模擬。ARMIE是ARM公司為了方便相關從業人員進行SVE指令的研究而推出的一個模擬器，ARMIE可以通過選項“-msve-vector-bits”指定模擬的向量寄存器寬度，支持從128位～2 048位共16種不同的向量寄存器寬度。受模擬器運行速度的限制，對測試用例完整計算過程的模擬難以在有限時間內完成，因此后續實驗針對每一個測試用例選取了計算過程中的部分時間迭代開展模擬。

(2)軟硬件平臺。ARMIE需要運行在支持ARM AArch64指令集的硬件平臺上。實驗所采用的是一個FT-2000+/64[15]的CPU，包含64個核，64 GB主存，實驗使用的編譯器為GCC-10.1.0，為了生成SVE指令，給編譯器提供的編譯選項為“-O2-march=armv8-a+sve”。

(3)SPH軟件平臺。本文選擇開源SPH項目DualSPHysics[16]作為本實驗的軟件平臺。DualSPHysics由西班牙維戈大學和英國曼徹斯特大學共同開發，采用C++語言編寫，適用于劇烈運動或重力場和液體的建模場景。本文對DualSPHysics項目中計算粒子間相互作用的部分進行了向量優化，粒子間相互作用的計算占據了DualSPHysics計算過程的大部分時間。對于粒子間相互作用的加速可以顯著提高整個軟件的性能。

(4)測試用例。實驗采用5個測試用例來驗證向量化的加速效果，對測試用例的簡要描述以及測試用例的規模(粒子數量)如表1所示。

對于上述5個測試用例，本文采用了2種不同的精度模式：單精度和雙精度，對比了從128位～2 048位16種不同向量寄存器寬度的情況下，粒子間相互作用計算部分，經過SVE向量化之后對比非向量化的加速比。

Table 1 Description of the testcases

單精度和雙精度2種不同模式下的加速比如圖7和圖8所示。可以看出雙精度模式下的加速比要高于單精度模式下的加速比。單精度模式和雙精度模式下的最高加速比來自于測試用例Pump，分別為3.71倍和4.92倍。加速效果與ARM針對其他Benchmark所給出的最高加速比存在差距的原因主要包括2個方面：一方面是因為粒子間相互作用的計算仍然存在串行計算部分，并不是所有的計算都可以通過向量化并行計算，比如由于粒子的位置發生了變化，需要更新粒子所在單元的信息；另外一方面，在算法2的SVE向量化算法中，第4行所遍歷的粒子均需在第7行計算其與粒子的距離，只有距離小于支持域半徑的粒子才會參與后續的向量運算，向量寄存器并沒有被填滿，這也是SVE的特點之一。

Figure 7 Speedup in single precision mode

Figure 8 Speedup in double precision mode

為了驗證SVE向量化算法的正確性，本文以測試用例DamBreakCubes為例，對比了非向量化算法的渲染結果與SVE向量化算法的粒子分布，結果如圖9和圖10所示。通過對比可以發現，經過SVE向量化加速之后的粒子分布與原始結果的渲染圖一致，表明本文所提出的SVE向量化方法是正確的。

Figure 9 Rendering effect of the testcase DamBreakCubes using serial algorithm

Figure 10 Particle distribution of the testcase CamBreakCubes using SVE algorithm

5 結束語

本文基于ARM SVE擴展指令集對光滑粒子流體動力學(SPH)開展了SIMD優化，對SPH計算過程中的粒子間相互作用計算進行了向量化并行，在單精度模式下取得了最高3.71倍的加速比，雙精度模式下取得了4.92倍的加速比。下一步將開展大規模SPH計算在SVE硬件平臺上的加速效果研究。