船舶三維聲彈性模擬軟件的并行優化策略*

呂小敬，劉 釗，蔣令聞，陳德訓，楊廣文1,

1.國家超級計算無錫中心，江蘇 無錫 214072

2.中國船舶科學研究中心，江蘇 無錫 214082

3.國家并行計算機工程技術研究中心，江蘇 無錫 214083

4.清華大學，北京 100084

1 引言

船舶三維聲彈性分析理論與方法研究彈性浮體與水介質的耦合振動及由此引起的聲輻射、聲散射和聲傳播問題。在此基礎上開發的船舶三維水彈性聲學分析軟件THAFTS-Acoustic（three-dimensional hydroelastic analysis of floating traveling structures）可以實現船內振動傳遞與船舶水中輻射聲場的統一計算與分析，具有良好的工程應用性[1-2]。

三維聲彈性力學的研究在改善船舶運動性能與安全性，控制船舶振動噪聲及提高水下隱身性能等一系列工程問題中有廣泛的應用需求和發展前景。1984年，Wu[3]建立了二維水彈性力學理論，將船體結構簡化為非均勻Euler 梁或Timoshenko 梁。Price 和Wu[4]將結構動力學理論與三維船舶運動勢流理論相結合，提出了廣義流固耦合邊界條件，開創性地發展了適用于波浪中任意三維可變形體承受內的三維水彈性力學理論。Du 等[5]發展了零航速三維脈動源Green函數快速計算方法，并建立了完善的三維航行船體線性水彈性力學頻域分析理論的數值計算方法。Zhou等[6]在三維水彈性理論及程序基礎上，發展了帶速航速、計及海面及海底邊界影響的船舶三維聲彈性理論，同時開發了一套較完整的可用于解決復雜船舶結構低中頻段聲彈性問題的數值模擬軟件。

三維聲彈性理論及軟件功能已日益完善，如何提高軟件計算能力，完善軟件計算復雜結構、復雜海洋水聲信道環境的功能，實現多工況、大任務計算的能力已迫在眉睫。而近年來，高性能計算蓬勃發展，如結合高性能計算理論及超級計算機的海量計算資源，對現有程序進行并行升級和優化，提高軟件的大規模高效率計算性能，成為了一個具有重要應用價值的研究課題[7-15]。

本文針對三維聲彈性理論算法特點，結合“神威·太湖之光”計算資源及體系架構，完成了三維聲彈性軟件THAFTS-Acoustic 廣義水動力系數計算模塊（THAFTS-Acoustic-hycof）核心算法的多級并行。使用某典型算例對多級并行優化策略進行測試和驗證，測試結果表明：完成多級并行后，核心算法代碼段獲得18 倍的加速效果，且程序具有良好的并行可擴展性，能夠支持THAFTS-Acoustic進行超大規模和更高精度的并行計算，為進一步改善我國船舶運動性能、安全性、震動噪聲以及水下隱身性能指出了一條可行的道路。

2 三維聲彈性算法及應用軟件簡介

2.1 三維聲彈性算法概述

頻域內船舶聲彈耦合動力學方程為：

其中，ω為頻率，MA為結構干模態廣義質量矩陣，CA為廣義阻尼矩陣，KA為廣義剛度矩陣，[AA]為干模態附連水質量矩陣，[BA]是干模態附連水質量阻尼矩陣，[CA]為廣義恢復力稀疏矩陣。

流場內的速度勢及基本運動方程滿足：

其中，?(x0,y0,z0,t)為相對于固定坐標系的流場中總的速度勢，為船舶恒速直線航行所誘導的相對于平衡坐標系的非均勻穩態流場速度勢，?O(x,y,z,t)、?D(x,y,z,t)、?R(x,y,z,t)分別為入射波速度勢、反射波速度勢及輻射波速度勢。

考慮到船體在無粘、無旋的可壓聲介質中，應滿足Bernoulli方程（Cauchy-Lagrange積分）：

結合各階模態的輻射波對應的流固耦合邊界條件：

采用頻率法求解時，r階廣義力F的表達式可展開為：

其中，Arj、Brj、Crj分別為流體作用在結構上的附連水質量矩陣、附連水阻尼矩陣及廣義恢復力系數矩陣，S為結構平均濕表面。

采用有限水深環境中聲波Green函數（7）求解流場內的速度勢及基本運行方程可得：

輻射波速度勢的表達式為：

采用Hess-Smith 等強度源方法，將船體濕表面（與流體接觸的表面）離散成N塊四邊形面元，則：

各階模態的輻射波速度勢與主坐標響應滿足：

根據模態疊加原理，由式（1）～式（10）求出各階干模態主坐標響應qr(r=1,2,…,m)后，帶入船舶結構振動方程，即得結構振動響應及流場內輻射聲性能。

三維聲彈性計算涉及多場耦合、多物理量、多核心段，單一的并行模式根本無法滿足所有計算熱點的高效并行，因此需要針對該軟件算法特點，構建多層次、多類型異構并行模型，支持數據并行與任務并行結合的混合并行模式，擴大程序并行度，保證各并行層次上的負載平衡，同時結合神威太湖之光超級計算機體系架構，研究恰當的多級混合并行算法實現，充分發揮眾核處理器的超高計算性能。

2.2 三維聲彈性軟件流程

三維聲彈性軟件包含三個模塊：flxbd、hycof、hyelas。flxbd模塊對輸入的數據進行預處理，生成廣義水動力系數計算模塊hycof所需數據；hycof模塊通過計算格林函數及其偏導數，計算源強和速度勢，得到水動力系數等參數；hyelas模塊根據水動力參數求解廣義流固耦合動力學方程，生成后處理所需數據。其中廣義水動力系數計算模塊hycof和hyelas計算量較大，目前僅實現了一維濕面元并行，程序并行效率低，無法滿足濕面元/模態平方依賴或者更高依賴的函數，且超過64進程時，會出現倒加速情況。本文針對hycof 模塊內不同計算過程及各計算過程數據依賴關系，考慮計算過程復雜度，考慮多級并行時各計算過程銜接及過程間數據銜接，實現了廣義水動力系數計算模塊多級異構并行，提高程序并行度及可擴展性能。THAFTS-Acoustic軟件流程圖如圖1所示。

Fig.1 Flow chart of THAFTS-Acoustic software圖1 THAFTS-Acoustic軟件流程圖

3 神威·太湖之光超級計算機系統簡介

3.1 SW26010處理器

神威·太湖之光超級計算機全機由40 960 塊SW26010 異構眾核處理器、20 480 塊計算板節點組成，共有10 649 600 個計算核心，系統峰值性能為125.4 PFlops，已連續4次蟬聯TOP500榜首。SW26010處理器架構如圖2所示。

SW26010異構眾核處理器包含4個核組，每個核組包含1 個主核和64 個從核，核組內提供8 GB 本地內存。從核局部（local data memory，LDM）存儲空間大小為64 KB。從核訪問LDM 速度較快，因此眾核優化的關鍵是減少從核訪主存次數，提高LDM 利用率。

3.2 “神威·太湖之光”并行方法

“神威·太湖之光”計算機系統語言環境包括基礎語言系統、并行編程語言接口、用戶使用環境及基礎編程環境，支持消息并行模型、共享并行編程模型、加速并行編程模型，同時支持4種異構并行方式：主從加速并行、主從協同并行、主從異步并行、主從動態并行。異構并行方式如圖3所示。

Fig.2 Architecture of SW26010圖2 SW26010處理器架構

Fig.3 Heterogeneous parallel method on Sunway TaihuLight圖3 “神威·太湖之光”計算機系統異構并行方法

主從加速并行中，計算核心通過加速線程庫加載核心段到從核上完成加速計算，主核等待加速任務結束后完成通信、IO（input/output）和部分代碼計算；主從協同并行時，主核作為一個核心，與從核一起完成核心任務計算；主從異步并行時，主核完成計算、通信、IO操作，從核完成核心加速，可實現計算通信隱藏、計算IO 隱藏，優化效果明顯；主從動態并行主要應用于從核計算任務時間不固定或者某些任務并行的程序。

4 三維聲彈性軟件多級異構并行算法實現及優化

本文基于異構并行模型，根據船舶三維聲彈性算法特點，以不同計算階段計算密度為特征，以變量依賴關系分析為基礎，在以數據預處理和基本變化為主的前后端計算階段采用MPI（message passing interface）消息傳遞并行；在核心計算階段，研究多級異構主從異步算法，增大程序并行度，采用分塊策略和通信計算隱藏策略，實現不同并行層次上的負載平衡。多級并行方案如圖4所示。

Fig.4 Multi-level parallel schematic of THAFTS-Acoustic-hycof圖4 三維聲彈性軟件多級并行圖

4.1 一維濕面元級并行算法實現

假定船體離散后總濕面元數為IXX，總進程數為numprocs，一維數據并行時數據并行方式如圖5所示。

Fig.5 One-dimensional parallel strategy of wet surface element圖5 一維濕面元級并行策略

采用一維濕面元并行操作較為簡單，對程序的修改最少，進程間通信次數和通信總的數據量可以達到最小，便于實現多級并行中各過程間數據銜接及多工況程序耦合。

4.2 二維濕面元并行及模態并行算法實現

一維濕面元并行可以較為簡單地實現計算量線性依賴濕面元數的函數，然而對于平方依賴或更高依賴濕面元的函數效果卻不是很好，如計算格林根數偏導VIN(IXX,IX)及求解源強SV(IXX,MODE)，隨著濕面元數及求解模態數增加，程序并行效率降低。以由格林函數偏導求解源強為例，同一列VIN分布在不同的進程中，主進程需要通信收集列主元行號，通信完成列主元行與當前行交換，完成處理后，仍需要與其他進程通信。當并行規模較小時，這種并行方案有一定并行加速效果，但是隨著進程數增多，通信量急劇增大，加速效果越來越差，甚至出現倒加速。因此需要設計更為合理的并行模式。

本文采用二維濕面元/模態分塊并行解決上述問題。假設行進程數為NPROW，列進程數為NPCOL，行進程分得塊大小為NNB，列進程分塊大小為NNM,各進程分得的行數為NNP，列數為NNQ，MYROW為行進程號，MYCOL為列進程號，MYID為當前進程組號，數據劃分方式如圖6所示。

二維分塊并行后，可以對進程進行分組操作，同時某進程完成NNB行×NNM列計算后，在與其他進程進行塊通信同時，可計算下一個NNB行×NNM列塊，便于實現通信和計算隱藏。

對于某些計算量與濕面元數/模態數成平方依賴的函數，采用二維濕面元并行程序擴展性更好，同時便于利用更大的進程數。

4.3 基于athread的眾核并行算法實現及優化

SW26010處理器上每個從核配備了用戶可控的64 KB局部數據高速緩存（LDM），支持gld/gst直接離散訪問主存及DMA批量數據訪問主存，并將數據放置在LDM 中。從核訪問主存訪問效率很低，需要數百個時鐘周期，而從核訪問LDM 數據僅需要數個時鐘周期。因此，眾核并行性能提升的關鍵是充分利用從核LDM，減少訪存次數，降低訪存開銷及通信開銷。

4.3.1 DMA通信數據的合并

Fig.6 Two-dimensional block parallel strategy of wet surface element/mode圖6 二維濕面元/模態分塊并行策略

減少通信和訪存開銷的最直接方法就是減少通信次數，減小需要通信的數據和增大單次數據通信的數據長度。因此可以將需要通信的數據進行合并和計算，減少通信次數，提高通信帶寬利用率。

三維聲彈性求解速度勢時，觀察到格林函數數組GRNN(IXX,5,IX)與AR_A(IXX)數組，按照LDM大小分塊拷貝到從核時，需要6次通信。將格林函數數組GRNN與AR_A數組合并為一個數組后，通信次數可降低為1，且可增大單次數據拷貝量，既減少了訪存，又增大通信帶寬利用率。

4.3.2 循環分塊及循環合并

異構眾核并行及優化中，需充分利用從核局存訪存性能，將數據盡量多地放到LDM 中，而SW26010 芯片的從核LDM 僅有64 KB，實際課題中無法將內層循環相關物理量一次性拷貝到LDM 中，訪主存則加速性能較低。本文根據從核LDM 大小，選擇合適的分塊并行方式，提高從核訪存性能。

實際應用課題中，核心段一般包含多重循環，如果僅按照最外層循環做任務劃分，易造成從核間負載不均衡。本文采用外層循環合并方式，增大程序并行度，提高程序的并行性能。循環分塊及循環合并實現如圖7所示。

Fig.7 Loop tile and loop collapse multicore parallelization scheme圖7 循環合并及循環分裂眾核并行方案

4.3.3 通信計算隱藏

對于眾核加速計算程序，提高DMA 帶寬，降低從核DMA 通信時間是提升性能的另一關鍵策略。本文采用計算和通信隱藏方案，具體實現過程如圖8所示。

Fig.8 Over hiding parallel scheme of communication and computation圖8 通信計算隱藏并行方案

4.3.4 向量化

SW26010處理器具備256位向量寄存器，可一次處理4 次浮點計算，8 次整型計算。因此一條SIMD（single instruction multiple data）指令相當于一個小的循環，可以減少指令數及由循環引起的控制相關，充分利用SIMD擴展結構提高性能。

SIMD編程在標準C的基礎上擴展得到了6種標準數據類型intv8、unitv8、int256、uint256、floatv4 及doublev4，標準變量需要通過顯式SIMD 內部函數調用、擴充的數據類型等實現SIMD的功能。

hycof 核心段內主要為復矩陣運算，計算兩個復數乘c=c+a*b時，假設初始c=(c1+c2i)，a=(a1+a2i),b=(b1+b2i),則c=(c1+a1*b1-a2*b2)+(c2+a1*b2+a2*b1)i。然而神威太湖之光系統沒有與復數類型匹配的擴展數據類型，必須將復數矩陣轉變為標準擴展數據類型。復數分實部虛部，一個complex*16 相當于兩個real*8。本文采用如下復數向量化方案，如圖9 所示：兩個復數可擴展為一個doublev4 擴展數據類型，之后通過SIMD 內部函數調用，可分離數據，得到復數計算實部虛部，再通過向量間運算，實現復數的運算操作，提高程序運算性能。

Fig.9 SIMD parallel scheme of three-dimensional acoustic elastic complex matrix圖9 三維聲彈性復矩陣SIMD并行方案

4.4 I/O優化

由于I/O 性能的增長速度跟不上系統本身處理能力的發展，I/O 性能成為高性能并行計算的主要性能瓶頸。THAFTS-Acoustic 軟件中采用文件保存臨時計算數組，當數據的讀寫量不大時，能夠滿足需求，但是隨著數據規模的增大，對I/O 提出了更高的要求。本文采用以下三個優化策略：（1）針對內存占較大的臨時數組，程序中的格林函數GRNN 采用全局變量在各過程間傳輸數據，刪除變量讀寫I/O；（2）數值模式中所有通信域內的進程均參與數據文件的輸出，各進程將其負責的數據寫入獨立文件，提高I/O 并發度；（3）進程分組收集數據，將各進程組內數據合并，減少寫文件I/O次數，增大單次I/O數據量。

5 并行性能測試及算例驗證

5.1 軟件并行性能測試

本文采用濕面元數5 678，模態數1 500，計算頻率數為3的左舷半球殼的算例。

5.1.1 多級異構并行加速性能

本文基于神威太湖之光系統，采用64進程及256進程兩個并行規模，對二維濕面元/模塊并行程序及多級異構并行程序相較一維濕面元并行程序的加速性能進行了初步的測試，如圖10所示。

64 進程測試時，二維濕面元/模態并行程序相較原始一維濕面元并行程序整體加速2.37 倍，多級異構并行程序相較原始一維濕面元并行程序整體加速5.54 倍。隨著進程數的增大，采用二維濕面元/模態并行及多級異構并行的優勢更為顯著。

本文采用4～512 進程測試多級異構并行程序的MPI 擴展性能，測試結果如圖11 所示。多級異構并行規模小于256 進程時，具有較為理想的加速效率，但是隨著數據規模的擴大，并行效率降低。這是由于算例濕面元數是5 678，當進程數大于256時，進程內數據量減少，計算開銷與通信、IO 開銷比明顯降低，并行效率下降。

Fig.10 Multi-level parallel acceleration performance圖10 多級并行加速性能

Fig.11 Strong scalability of multi-level heterogeneous parallelization圖11 多級異構并行強可擴展性測試

5.1.2 眾核加速性能測試

本文采用濕面元數5 678，模態數1 500，計算頻率數3 的左舷半球殼的算例。基于神威太湖之光系統，采用256進程測試眾核并行加速性能。測試的對象包括源強求解及速度勢求解兩個過程。程序眾核加速性能如表1所示。

Table 1 Multicore parallel acceleration ratio of kernel functions in THAFTS-Acoustic-hycof表1 THAFTS-Acoustic-hycof核心函數眾核并行加速比

采用向量化、通信與計算隱藏、循環合并及循環分裂等優化方法完成眾核并行后，程序求解速度勢加速17.7 倍。源強求解核心函數以塊為單位并行，求解中存在大量的進程間通信，在計算完成后需要將數據塊使用進程組通信發送給其他進程，而通信過程無法使用眾核進行加速。在通信成為瓶頸的情況下，源強求解函數仍然獲得了6.6 倍的整體加速效果。

5.2 軟件正確性驗證

本文采用濕面元數5 678，模態數1 500，計算頻率數3的左舷半球殼的算例，驗證多級異構并行程序正確性。3個頻率點（每個頻率點計算2次，共計算6次）的物理量附連水質量、附連水阻尼相較原始程序的協方差如圖12所示。

由圖12可知：6次計算的附連水質量協方差均在E-16 量級，附連水阻尼的協方差均在E-14 量級，誤差范圍合理。多級異構并行結果與軟件初始版本計算結果一致。

本文采用經典的受徑向單位集中力作用下舷間充水雙層彈性球殼結構聲輻射的算例模型進行軟件正確性驗證。計算模型如圖13 所示：該計算模型存在解析解及THAFTS 軟件模擬解，具體可參考文獻[6]。

采用零航速，無界流場，內球殼半徑0.5 m，內球殼壁厚1 mm，外球殼半徑0.65 m，外球殼壁厚0.3 mm，內外球殼體密度7 800 kg/m3，楊氏模量2.1×1 011 N/m2，泊松比0.3，干模態阻尼比0.01，內外場流體密度1 025 kg/m3，內外場流體聲速1 500 m/s。內球殼濕面元網格邊長約0.05 m，外球殼濕面元網格邊長約0.065 m。計算4個場點的輻射聲壓。軟件求解和解析計算結果比對情況如圖14所示。

Fig.12 Correctness verification圖12 正確性驗證

Fig.13 Double layer concentric sphere and coordinate system圖13 雙層同心球和坐標系

測試結果表明，多級并行后，三維聲彈性計算程序是正確的，可應用到實際工程領域。

Fig.14 Comparison of calculation results圖14 計算結果對比

6 總結與展望

三維聲彈性力學的研究可廣泛應用在改善船舶運動性能與安全性的應用領域。在本文中，完成了三維聲彈性核心求解模塊THAFTS-Acoustic-hycof的多級異構并行工作。主要包括以下幾點：根據程序算法復雜度，完成了THAFTS-Acoustic-hycof 的一級濕面元并行、二級濕面元/模態并行，設計了三維聲彈性源強及速度勢求解過程的眾核并行及優化，包括DMA 通信數據合并、循環分裂及循環合并、通信計算隱藏及SIMD等。測試結果表明，三維聲彈性多級并行程序具有良好的并行加速比和眾核加速性能，能夠有效地發揮SW26010 國產眾核處理器和神威?太湖之光超級計算機的強大計算能力，從而大大縮短實際工程應用的項目周期，加速改善我國船舶運動性能、安全性、震動噪聲等關鍵性能，提高船型整體設計研發效率。

未來將在以下方面開展工作：

（1）THAFTS-Acoustic-hycof 軟件的頻率并行：目前測試階段計算頻率個數有限，但潛在頻率數比較大，因此如何實現頻率并行與濕面元、模態并行，眾核并行結合是進一步提高程序并行效率的關鍵。

（2）斷點恢復及容錯功能：完善THAFTS-Acoustichycof在大規模異構并行計算機系統硬件故障容錯和重啟動功能，保證計算結果的正確性和可靠性。

（3）THAFTS-Acoustic 軟件的全過程并行及優化：開發THAFTS-Acoustic-flxbd 前處理模塊并行算法及程序，開發THAFTS-Acoustic-hyelas后處理模塊的多級并行及優化。