三維水彈性軟件并行化的初步實現(xiàn)

葉永林，田 超，張彥斌

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

0 引 言

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

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

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

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

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

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

1.1 基本理論

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

1.2 程序框架與計算熱點

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

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

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

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

第三個程序HYELAS求解各階模態(tài)的主坐標響應(yīng)，其計算量不太大，所用時間較短。Hyelas函數(shù)運行時間統(tǒng)計如圖2所示。

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

從上面統(tǒng)計結(jié)果看，由INFLCO函數(shù)及調(diào)用子函數(shù)計算比重很大，而hyelas從目前結(jié)果看87.5%為系統(tǒng)調(diào)用，二者GAUSS1消元所占比重都有所降低。

2 并行計算方法簡介

2.1 并行計算原理

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

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

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

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

2.2.1 潛在計算量分析

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

分析程序，可看到hycof程序有多個并行，即數(shù)據(jù)計算循環(huán)較多，其中計算量最大的為入射波頻率、結(jié)構(gòu)濕面元并行。首先應(yīng)進行的、且最主要的并行級為濕面元的并行。而計算頻率個數(shù)目前規(guī)模有限，但潛在頻率數(shù)會比較大，因此本次并行中也需要頻率級并行實現(xiàn)。

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

對于hycof程序子函數(shù)FLXHYFP，還存在平方依賴于結(jié)構(gòu)彈性模態(tài)自由度。

2.2.2 軟件并行整體框架設(shè)想

實現(xiàn)hycof程序在濕面元并行框架方案如下：

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

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

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

3 程序并行實現(xiàn)

3.1 并行化流程

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

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

3.2 hycof分段并行化實施方案

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

（a）INCIDWP

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

（b）INFLCO

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

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

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

（c）SOLVQ1P/SOLVQ1及 SOLVPDP/SOLVPD

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

調(diào)用子函數(shù)方面，F(xiàn)04ADF（GAUSS1）需要后面獨立闡述。PHIDBC函數(shù)中有單重循環(huán)，因此循環(huán)和相關(guān)數(shù)組變量在濕面元維切分。BOUND函數(shù)類似。READISK函數(shù)經(jīng)過整理和優(yōu)化后，可讀性和性能都獲得提升。

（d）GAUSS1

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

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

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

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

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

（e）FLXHYFP及 FLXHYF

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

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

（f）DIFORCEP 及 DIFORCE

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

3.3 頻率級并行

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

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

4 算例效果

4.1 正確性驗證

4.1.1 常規(guī)船舶

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

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

使用原水彈性力學分析程序及MPI并行程序?qū)ι鲜龃霸趩挝徊ǚ?guī)則波下的船體主坐標響應(yīng)進行計算，其測試結(jié)果如圖5所示。其中，mpi為并行化程序計算結(jié)果。

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

圖5 主坐標響應(yīng)計算結(jié)果對比Fig.5 The parallel computing procedure of hydroelasticity code for VLFS’response

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

圖7 主坐標響應(yīng)計算結(jié)果[13]Fig.7 Principle response of VLFS in wave

4.1.2 超大型浮體

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

結(jié)果表明，并行化三維水彈性分析程序?qū)Τ笮透◇w的響應(yīng)分析結(jié)果正確，程序適用。

4.2 快速性效果

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

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

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

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

5 結(jié)論與展望

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

[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船結(jié)構(gòu)振動與噪聲分析[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]崔維成.超大型海洋浮式結(jié)構(gòu)物水彈性響應(yīng)預報的研究現(xiàn)狀和發(fā)展方向[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]路 振,倪歆韻.復雜形狀超大型浮式結(jié)構(gòu)物單模塊水彈性響應(yīng)研究（Ⅱ）[R].無錫:中國船舶科學研究中心科技報告,2014.Lu Zheng,Ni Xinyun.Hydroelastic response of single module of VLFS with complex shape(Ⅱ)[R].Wuxi:Report of CSSRC,2014.

[7]田 超,丁 軍,楊 鵬.島礁波浪環(huán)境下浮式結(jié)構(gòu)物的動響應(yīng)預報[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]金建海,葉永林,田 超.三維水彈性集成軟件的設(shè)計與實現(xiàn)[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].北京:電子工業(yè)出版社,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]田 超,陳映秋.基于三維水彈性力學的船體振動響應(yīng)計算[C]//北京力學會第16次學術(shù)年會,北京,2009.Tian Chao,Chen Yingqiu.Structural vibration response by hydroelasticity method[C].16th Conference of Beijing Mechanics Institute,Beijing,2009.

[12]超大型浮體結(jié)構(gòu)設(shè)計圖紙[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.