三維水彈性軟件并行化的初步實現

葉永林，田 超，張彥斌

（1中國船舶科學研究中心，江蘇 無錫 214082；2無錫市江南計算技術研究所，江蘇 無錫 214082）

0 引 言

海洋浮體在風浪的作用下會產生剛體運動，同時還發生彈性變形，彈性變形對周圍流場會有影響，從而也會改變波浪作用在結構上的力。水彈性理論概念最早于二十世紀五十年代從航空領域引入海洋工程領域，隨著耐波性理論中二維切片理論的引入，分析彈性船體對稱響應的二維水彈性理論和分析彈性船體反對稱響應的二維水彈性理論才出現[1]。上世紀八十年代中期，吳有生[2]開創性地將三維適航性理論與三維結構動力學理論相結合，提出了廣義流固界面條件，發展了一個適用于分析波浪中任意三維可變形體承受內、外激勵時動響應性能的三維水彈性理論。

船舶水彈性力學從動力學的觀點出發，將柔性船體與周圍流場作為一個相互作用的整體系統進行描述和分析，并在此基礎上預報船體結構的動力學特性，諸如：船體運動、波浪載荷、結構變形、應力應變以及疲勞性能等，以更合理地預報和評估船舶及海洋工程結構的運動特性、可靠性和安全性。現有的三維水彈性分析方法已經應用于上述諸多領域，取得了很好的效果[3-4]。

超大型海上浮式結構物的結構形式極其復雜，三維尺度龐大，在風浪流環境長期作用下和風暴氣象條件下可能產生的結構破壞模式眾多，遠遠超出了船舶和一般海洋工程結構物的范疇。對于尺度以數百到數千米計的超大型海上浮式結構物來說，由于動力載荷所產生的水彈性運動和結構變形不能被忽略，因此對超大型海上浮式結構物來說，采用水彈性理論計算并考慮各模塊間的水動力影響以及流固耦合問題是必須的[5]。

而近島礁超大型浮式結構物不僅形狀復雜、尺度巨大，同時結構物所處環境為非均勻海底淺水環境，因此，需要用到大量的面元數去離散浮式結構和非均勻海底。較大的面元數意味著所獲得的影響系數矩陣為一個階數巨大的非對稱復系數稠密矩陣，此外，較多的波浪頻率和模態數都導致了龐大的計算量，分析過程十分緩慢，急需開發并行計算方法，提高程序的運行效率[6-7]。

本文的工作主要介紹在中國船舶科學研究中心現有的三維線性頻域水彈性分析程序 THAFTS基礎上，通過對程序的計算熱點分析，形成新的程序框架和數據傳輸方法等，形成新的具有大型數據計算能力的快速并行化程序。并通過常規船舶與超大型浮式結構物的算例，驗證了并行化程序的正確性，測試了大型數據計算能力和速度。

1 三維線性水彈性分析簡介

1.1 基本理論

假定浮體結構周圍為理想、不可壓縮流體，波浪幅值為小量，可得到浮體的線性水彈性力學運動方程為[2]：

1.2 程序框架與計算熱點

現有的三維水彈性系列程序都是以吳有生院士的三維線性頻域水彈性力學程序為基礎發展起來的。目前三維線性頻域水彈性程序包THAFTS由三個模塊程序組成，程序語言為Fortran[8]。

圖1 Hycof程序計算熱點Fig.1 Calculation hot of Hycof program

其中第一個程序FLXBD計算彈性結構濕表面幾何特性以及各階模態的濕表面邊界條件和廣義流體恢復力系數；該部分計算量不大，所用時間較短。

第二個程序HYCOF占了程序的絕大部分計算時間，該項程序計算彈性浮體結構周圍流場速度勢，計算結構在單位幅值規則波中剛體運動和彈性變形模態的廣義水動力系數和廣義波浪激勵力。Hycof計算20個頻率運行時間比超過1%的所有函數如圖1所示。

第三個程序HYELAS求解各階模態的主坐標響應，其計算量不太大，所用時間較短。Hyelas函數運行時間統計如圖2所示。

圖2 Hyelas計算熱點Fig.2 Calculation hot of Hyelas program

從上面統計結果看，由INFLCO函數及調用子函數計算比重很大，而hyelas從目前結果看87.5%為系統調用，二者GAUSS1消元所占比重都有所降低。

2 并行計算方法簡介

2.1 并行計算原理

20世紀70年代，向量機的出現標志著現代超級計算機的快速發展，它使得計算機性能大大提升。80年代，考慮到并行可擴展性，基于分布式存儲（即各個處理器都有自己獨立的存儲設備）的并行計算機成為研究熱點。21世紀以來，超大規模并行機（MPP）以更優越的性能進入工程與科研領域，在航空航天、生物化學、地球物理等領域取得了大量的應用[9]。

并行編程環境方面，即并行計算機系統編程的特定抽象，是指為構建并行程序提供的基本工具、語言功能和應用編程接口。隨著并行計算機的發展，到20世紀90年代中期，出現了上百種并行編程環境。這些并行環境都基于如下并行體系結構：共享存儲器、具有消息傳遞的分布式存儲器、或兩種體系混合[10]。目前并行編程社區匯聚于兩個主要的并行編程模型：OpenMP和MPI。

OpenMP是一些語言擴展的集合，并被實現為編譯器指令，目前支持C、C++、Fortran，主要用于共享存儲系統上的單進程并行性的開發。MPI是一些例程的集合，這些例程提供進程管理、消息傳遞和某些進程相關的通信操作。MPI支持多種語言，主要用于消息傳遞的分布式存儲器并行計算機系統。本文采用了MPI環境進行三維水彈性程序的并行化改造。

2.2 三維水彈性程序并行可行性

2.2.1 潛在計算量分析

通過軟件具體分析和實際測試，發現目前三維水彈性程序THATFS的3個程序中，flxbd獨立性較強且計算量不大，暫時不需并行，計算量最大的是hycof，而hyelas有一定的潛在計算量。

分析程序，可看到hycof程序有多個并行，即數據計算循環較多，其中計算量最大的為入射波頻率、結構濕面元并行。首先應進行的、且最主要的并行級為濕面元的并行。而計算頻率個數目前規模有限，但潛在頻率數會比較大，因此本次并行中也需要頻率級并行實現。

程序計算量方面，對hycof內部運行部分而言，所有函數的計算量基本上依賴于濕面元自由度，大多數函數線性依賴于濕面元數，其余是平方依賴或更高。因此，當浮體濕面元增加時，整個程序的計算量將非線性增加。其中可以預見，以GAUSS1函數導致計算量增加最為明顯。

對于hycof程序子函數FLXHYFP，還存在平方依賴于結構彈性模態自由度。

2.2.2 軟件并行整體框架設想

實現hycof程序在濕面元并行框架方案如下：

設定參數IXX的值為濕面元數NP。假定有numprocs個MPI進程，設定IX=(NP+numprocs-1)/numprocs)，則前面（numprocs-1）個MPI進程切分的NP濕面元數NP_p為IX，而最后一個進程切分相應的濕面元數NP_p為NP-myid*IX。程序中循環變量范圍從1～NP切分為1～NP_p，相應的計算量減少。

在內存使用方面，除了個別變量外的大部分變量都在濕面元維度上進行了數據切分，尤其是原來與濕面元成平方比的數組，如VIN（IXX，IXX），在其中一個濕面元維度進行切分VIN（IX，IXX）。切分后的變量占用內存大小反比于運行的MPI規模，即計算所使用MPI進程規模越大，內存占用越小。軟件整體內存，也隨MPI規模增加而減小并最終趨向于一個小的固定值。

另外，hycof上程序中有大量的I/O和文件操作。在并行化時除了個別占用內存較大的數組之外，其余相對較小的數據都沒有采用0進程讀入再通過MPI通訊分發，而是需要每個進程讀入。這樣做基于如下考慮：并行規模不大時，程序性能表現更好，同時軟件編程實現和可讀性都更好。數據輸出部分在并行化時原則上保持了輸出文件不變，即文件輸出由其中一個主進程—即0進程來實現文件輸出。這樣做方便計算結果后處理。

3 程序并行實現

3.1 并行化流程

開展近島礁超大型浮式結構物的水彈性響應分析時，針對多級并行編程模型的并行算法設計與高效實現，通過濕面元、頻率等各級向量層的多級并行，并在此基礎上實現進一步的線程級并行，保證大規模并行重點應用的并行效率和整體效率具體工作流程為圖3所示。

圖3 超大型浮體水彈性響應并行化工作流程Fig.3 The parallel computing procedure of hydroelasticity code for VLFS’response

3.2 hycof分段并行化實施方案

因三維水彈性程序并行重點在于Hycof程序，因此本節就hycof各個具體程序段按照調用次序進行闡述：

（a）INCIDWP

本程序僅有一個循環，且結構較為簡單，輸入數組變量為PC，而每個進程都讀入了PC，且該數組變量未切分，輸出數組變量TIWP(IX)、EKX(IX)和EKY(IX)，都是切分后的數組變量。只需要注意循環上下界即可。另外，程序中多次調用INCIDWP函數，改為僅主函數調用一次即可。

（b）INFLCO

INFLCO函數較為復雜，其計算量也較大。圖1所述的熱點函數表中，有多個函數屬于INFLCO函數調用的子函數。

該函數有一個兩重循環，其計算量正比于濕面元數的平方。主要輸入變量PN、PC、AR等，輸出變量為GRNN、GRSS、VIN和VINS。對于這里存在兩重循環，在外層I循環進行濕面元切分。由于INFLCO調用的子函數中仍然有PN、PC、AR等變量，因此注意一些子函數調用形參中的I需要替換為新的 II=IX*myid+I。 VIN（IXX，IXX）也替換為 VIN（IX，IXX），依次類推。

特別說明的是，基于性能提升的考慮，將本程序中格林函數GRNN、GRSS進行提升維度處理，即將GRNN(IXX,4)替換為GRNN(IXX,4,IX)。提升維度后，GRNN和GRSS變量將常駐內存中，不必再寫回文件，也不必在后續使用時再次從文件中讀入，即去掉格林函數相關的I/O操作。

（c）SOLVQ1P/SOLVQ1及 SOLVPDP/SOLVPD

SOLVQ1/SOLVQ1P被SOLVE/SOLVEP調用，兩者結構基本一致。其主要輸出變量BV、QQS、QQA在濕面元維切分，這里GRNN、GRSS和AR連續維（即第一維）是全濕面元。程序中由ISOURCE參數控制兩條運行分支路徑，當ISOURCE為2時，有兩重循環計算VIN和BV，這里計算去掉了原來的if分支，程序更簡潔，循環在最外層IU進行切分。而ISOURE為1時則直接讀取VIN。在經過GAUSS求解后，ISOURCE為2分支較簡單，ISOURCE為1分支，存在三重循環，仍在最外層IU進行切分。

調用子函數方面，F04ADF（GAUSS1）需要后面獨立闡述。PHIDBC函數中有單重循環，因此循環和相關數組變量在濕面元維切分。BOUND函數類似。READISK函數經過整理和優化后，可讀性和性能都獲得提升。

（d）GAUSS1

軟件中GAUSS1比重不大，暫時僅考慮和其它子函數對應數組良好銜接性：

對于待求解行列矩陣A按照行切分，對于矩陣大小為N，使用有效進程數N_lastid，最后一個有效進程分得行數為N_p=N-myid*IX，其余N_p=IX。

出于負載平衡的考慮，列主元的確定由0進程負責，因為0進程的行消元過程最早完成。同一列的數據分在不同的進程，需要通信收集至0進程；獲得最大列主元行I0后廣播或通信出去，再通過通信實現I0行和當前行K進行交換；當前行進行對角元素歸一化的處理后再發給后面的進程；后面進程并行進行消元處理。

回代過程從N行開始，求得B（N）的解后，通過通信過程將B（N）發給前面的進程，其它進程并行回代。

回代過程完成后，通信收集所有解SV再廣播其它進程，后面其它函數需要整個濕面元場求和。

（e）FLXHYFP及 FLXHYF

本函數除調用的子函數以外的程序段不存在濕面元維的循環，只將相關的變量切分即可。在其調用的子函數BOUND1需要對循環切分。在AB函數中由于存在全濕面元空間求和情形，因此要針對濕面元維循環切分，再將待求和變量進行MPI_REDUCE操作。PRESSURE函數類似，獲得全濕面元空間變量需要MPI_GATHER并由0進程完成輸出。

FLXHYF的關于HYCOF.SOL及其它I/O部分由0進程負責。特別需要說明，該I/O部分運行開銷在模態較大時也會明顯增加，成為MPI加速的瓶頸之一，未來振動模態數增加時改為二進制格式輸出。目前振動模態數較少，為了保證不改變輸出文件格式和內容，暫時不優化該部分。

（f）DIFORCEP 及 DIFORCE

DIFORCEP函數調用SOLVPDP函數計算QDP、QQS、QDN和QQA變量后，再計算HS，最終對全濕面元維度求和獲得輸出的FH、FII和FDD變量。因此除了對所有變量和濕面元循環體在濕面元維度進行切分外，還要將相關變量在局部濕面元空間求和，然后進行MPI_REDUCE操作。I/O部分仍然由0進程負責。

3.3 頻率級并行

目前規模下，濕面元級的并行規模增加時，通訊比重會上升。在計算資源足夠的情形下，可以進一步考慮入射波頻率級并行，降低濕面元級并行規模。頻率級并行方案如下：將給定的numprocs個進程，全局的進程編號myrank，均分成numgroup個進程組，組編號為groupid，每個組有groupprocs個進程，組內每個進程編號為myid。每個組負責處理若干獨立連續頻率點的計算，進程組間數據無相關性。

頻率級并行的進程分組主要通過調用MPI_GROUP_INIT子程序中的MPI_COMM_SPLIT函數來實現。頻率級并行實現過程中對I/O部分的處理是關鍵。Hycof程序運行過程中的臨時文件，即TAPE為開頭的文件，以TAPE*myrank標號命名。最終的輸出文件如HYCOF.SOL等必須在頻率循環外以HYCOF.SOL*groupid，并且此類文件由每個進程組中myid為0的進程負責輸出。計算完成后，可以通過命令來實現文件的拼接。

4 算例效果

4.1 正確性驗證

4.1.1 常規船舶

首先采用18萬噸散貨船，進行并行程序的正確性測試，相關數據如圖4所示[11]。

圖4 散貨船主尺度Fig.4 Main dimensions of the sample ship

使用原水彈性力學分析程序及MPI并行程序對上述船舶在單位波幅規則波下的船體主坐標響應進行計算，其測試結果如圖5所示。其中，mpi為并行化程序計算結果。

常規船舶的水彈性響應分析結果表明，并行化后的軟件計算結果與串行軟件解算結果一致，計算結果正確合理。

圖5 主坐標響應計算結果對比Fig.5 The parallel computing procedure of hydroelasticity code for VLFS’response

圖6 超大型浮體算例主尺度[12]Fig.6 Main dimensions of VLFS

圖7 主坐標響應計算結果[13]Fig.7 Principle response of VLFS in wave

4.1.2 超大型浮體

本節針對如圖6所示的超大型浮體單模塊，進行了并行化計算測試，計算結果為單位波幅規則波時迎浪下的結構主坐標響應，結果如圖7所示。

結果表明，并行化三維水彈性分析程序對超大型浮體的響應分析結果正確，程序適用。

4.2 快速性效果

對上述超大型浮體單模塊的水彈性分析時，改變環境變量，測試計算速度。結果表明，當結構濕面元數較小時，僅用少量CPU，即可取得較好的計算加速效果；此時CPU增加，計算加速效果不會同比增加，如下表中的計算工況2、3所示。當計算入射波頻率數較多時，采用入射波頻率并行優化，可以取得更好的計算加速效果，如計算工況4所示。

對于超大規模算例（如圖8所示的多模塊浮體），測試結構的濕面元數為419 200個，此面元數基本可用于模擬超大型浮體的復雜構型與地形。此時其計算量非線性增加，計算量超大，串行計算異常緩慢，無計算結果可比。在測試平臺上只測試濕面元級并行軟件，當14個CPU時，單個波浪頻率的計算耗時為23h51m16.556 s。實際工程計算時，濕面元數將大幅減小，此時可增加CPU數量，使工程分析時間在較短時間內完成。

圖8 3模塊超大型浮體Fig.8 VLFS with 3 modules

表1 快速并行化效果（面元數：721）Tab.1 Effect of parallel computing(panel No.:721)

5 結論與展望

本文在中國船舶科學研究中心現有的三維線性頻域水彈性分析程序THAFTS的基礎上，根據超大型浮體在近島礁環境下的響應分析需求，形成快速計算的三維水彈性并行程序，并針對常規水面船及超大型浮體，進行了程序的正確性與快速性測試，結果表明，該并行程序與原程序計算結果一致。計算速度隨著CPU的增加而增加，在計算數據巨大時，加速性能突出。下一步將根據超大型浮體計算分析需求，開展彈性模態級并行，并進行濕面元級、入射波頻率級等的二級并行優化，使三維水彈性力學分析程序更加快速、高效、穩定，該程序在超大型浮體的響應分析中將發揮重要作用。

[1]Bishop R E D,Price W G.Hydroelasticity of ships[M].London:Cambridge University Press,1979.

[2]Wu Y S．Hydroelasticity of floating bodies[D].Brunel University,1984.

[3]田 超,吳有生.船舶水彈性力學理論的研究進展[J].中國造船,2008,49(4):1-11.Tian Chan,Wu Y S.Review of research on the hydroelasticity of ship[J].Shipbuilding of China,2008,49(4):1-11.

[4]葉永林,吳有生,鄒明松.基于水彈性力學的SWATH船結構振動與噪聲分析[J].船舶力學,2013,17(4):430-438.Ye Yonglin,Wu Yousheng,Zou Mingsong.Research on application of hydroelasticity method to vibration response analysis of the SWATH ship hull[J].Journal of Ship Mechanics,2013,17(4):430-438.

[5]崔維成.超大型海洋浮式結構物水彈性響應預報的研究現狀和發展方向[J].船舶力學,2002,6(1):73-90.Cui Weicheng.Current status and future directions in predicting the hydroelastic response of Very Large Floating Structures[J].Journal of Ship Mechanics,2002,6(1):73-90.

[6]路 振,倪歆韻.復雜形狀超大型浮式結構物單模塊水彈性響應研究（Ⅱ）[R].無錫:中國船舶科學研究中心科技報告,2014.Lu Zheng,Ni Xinyun.Hydroelastic response of single module of VLFS with complex shape(Ⅱ)[R].Wuxi:Report of CSSRC,2014.

[7]田 超,丁 軍,楊 鵬.島礁波浪環境下浮式結構物的動響應預報[J].船舶力學（已錄用）,2014.Tian Chao,Ding Jun,Yang Peng.Prediction of dynamic response for floating structures under wave environment around islands and reefs[J].Journal of Ship Mechanics,2014.

[8]金建海,葉永林,田 超.三維水彈性集成軟件的設計與實現[J].船舶力學,2012,15(5):521-529.Jin Jianhai,Ye Yonglin,Tian Chao.Design and implementation for 3D hydroelasticity integrated software[J].Journal of Ship Mechanics,2012,15(5):521-529.

[9]莫則堯,譯.并行計算綜述[M].北京:電子工業出版社,2005.Translated by Mo Zeyao.Summarizations of parallel computation[M].Beijing:Electronic Industries Press,2005.

[10]熬富江.并行編程模式[M].北京:清華大學出版社,2008.Ao Fujiang.Mode of parallel program[M].Beijing:Tsinghua University Press,2008.

[11]田 超,陳映秋.基于三維水彈性力學的船體振動響應計算[C]//北京力學會第16次學術年會,北京,2009.Tian Chao,Chen Yingqiu.Structural vibration response by hydroelasticity method[C].16th Conference of Beijing Mechanics Institute,Beijing,2009.

[12]超大型浮體結構設計圖紙[R].無錫:中國船舶科學研究中心科技報告,2013.Blue Print of VLFS’Hull[R].Wuxi:Report of CSSRC,2013.

[13]葉永林.超大型浮式單模塊水彈性并行分析[R].無錫:中國船舶科學研究中心科技報告,2014.Ye Yonglin.Hydroelastically parallel computation on single module of VLFS[R].Wuxi:Report of CSSRC,2014.