基于Woodbury+OpenMP的結構非線性地震反應并行分析方法

余丁浩， 李 鋼

(大連理工大學 建設工程學部海岸和近海工程國家重點實驗室， 遼寧 大連 116024)

隨著工程結構逐漸向大型化和復雜化方向發展，基于傳統變剛度法的結構非線性地震反應計算平臺所需巨大計算耗時已越發難以滿足工程建設需求，盡管目前計算機技術已獲得了長足發展，但大型結構抗震分析所面臨的效率瓶頸依然是制約結構工程發展的重要因素，如何在保證高精度計算模擬的同時有效提升結構分析效率始終是工程界重點關注問題。

利用計算機的多線程計算能力或搭建集群化計算環境通過并行方式進行結構非線性地震反應分析是提升計算效率的重要途徑，從編程角度來看可將當前結構分析領域常用的并行計算模式分為兩類：共享內存式和分布內存式。共享內存式是指多個計算線程共用一個數據存儲空間，主要用于單機計算，如基于OpenMP的多核CPU(central processing unit)并行計算[1-2]，基于GPU的并行計算[3-11]等，具有維護和編程相對簡便的優點；分布內存式主要利用計算機集群實現高性能計算，不同計算設備之間通過必要的數據通訊進行協調，此類并行計算通常采用消息傳遞的編程模式(message passing interface，MPI)實現，雖然其理論上可以實現超大規模問題的計算分析，但編程方式復雜、設備維護代價高，且不同計算設備之間無法避免的數據通訊導致計算效率存在瓶頸。在實際應用中，上述兩種并行模式可以混合使用，以進一步降低分析耗時，例如何強等[12]通過對整個求解域進行剖分，將不同子域分配到不同計算設備上，并在單機上利用OpenMP進行線程級別并行加速，提出了基于MPI+OpenMP混合編程模式的并行計算方法。此外，為充分發揮不同并行計算模式的效能，眾多學者有針對性的提出了多種便于并行實現的結構分析方法，如單元接單元方法[13-16]、區域分解法[17-21]等，其中單元接單元方法不需總裝整體剛度矩陣和荷載向量，運算主要集中在單元層面，該方法通常與預處理共軛梯度法相結合，多用于共享內存式的并行計算；區域分解法采用分而治之的思想，將結構劃分為不同的計算區域，其關鍵在于建立能夠實現負載平衡的區域劃分方式和界面方程求解方法，多用于分布式計算集群。

結構在地震作用下的非線性變形通常集中于局部區域，如部分梁柱端部、墻底等部位，已有研究表明，充分利用結構的局部非線性特征有助于大幅提升非線性地震反應分析效率，當前已發展出多種高效局部非線性分析方法，如擬力法[22-26]、數值子結構法[27-30]、多尺度類方法[31-32]、快速非線性分析法[33-34]等，這些方法多通過保持結構剛度矩陣恒定或縮減問題規模的方式避免結構大規模剛度矩陣實時更新及由此導致的計算效率降低。在現有的局部非線性分析方法中，基于Woodbury公式的結構非線性分析方法[35-37](簡稱Woodbury方法)因具有較好的數學理論基礎和較廣泛的適用性，近年來得到了越來越多的關注，此類方法將任意時刻的結構切向剛度改寫為初始彈性剛度的低秩攝動形式，進而引入線性代數領域中的Woodbury公式求解結構切向位移反應，能夠避免了大規模整體切向剛度的迭代更新，僅需對代表局部非線性行為的小規模矩陣進行實時更新和分解，可顯著提升非線性問題計算效率。Deng等通過將結構切向剛度中局部時變元素單獨取出，建立了以Woodbury公式為高效求解器的虛擬力分析法。近期Li等[38-40]以擬力法中的變形分解思想為出發點，通過對局部區域內進入塑性的材料進行應變分解，在單元層面構造非線性應變場函數，提出了隔離非線性有限元法基本理論，該方法采用Woodbury公式進行控制方程求解，且具有適用性廣、通用性好等優點，能夠實現一般結構的高效非線性地震反應分析。隨后，眾多學者圍繞該方法的算法優化、計算應用等開展了系列研究[41-48]，進一步豐富了其理論體系。盡管Woodbury方法能夠從算法理論層面解決傳統非線性分析方法由于剛度矩陣時變導致的效率瓶頸，但現有研究均采用串行計算模式，并未充分利用計算硬件的并行加速能力，因而其高效性優勢并未得到充分發掘。

本文將OpenMP并行計算模式引入到Woodbury方法中，提出了一種基于Woodbury+OpenMP的高效結構非線性地震反應分析方法。首先采用隔離非線性法構造結構切向剛度方程的攝動展開表達式及相應Woodbury求解公式，隨后，分別針對Woodbury公式計算、單元狀態確定、局部非線性系數矩陣更新這3個非線性迭代求解過程中主要計算步驟建立了基于OpenMP的并行優化策略，實現了該方法的全過程并行加速，最后，通過一個高層結構算例驗證了本文方法的高效性和準確性。

1 Woodbury方法基本理論

本文采用隔離非線性法構造用于結構地震反應分析的Woodbury求解公式。使用圖1所示結構對該方法中控制方程的基本推演過程進行說明，設在某個計算步下單元i進入非線性，圖1給出了該單元中任一點材料的應力應變關系，其中σ和ε分別為當前的材料應力和應變，根據材料初始彈性模量Ee可將應變ε分解為彈性和非線性兩部分，其中彈性應變ε′等于假設材料保持初始彈性狀態時達到應力σ對應的應變(即ε′=σ/Ee)，非線性應變ε″等于總應變ε與彈性應變ε′的差(即ε″=ε-ε′)。進一步利用插值方法在單元層面構造材料非線性應變分布場函數，并結合虛功原理，即可建立增量形式的單元隔離非線性控制方程

(1)

圖1 隔離非線性法基本原理

(2)

在任一迭代計算步中，將式(2)中的非線性自由度進行凝聚可得攝動式剛度方程式(3)。

(3)

式中:n為結構位移自由度數；m為當前計算步下結構的非線性自由度數，代表了結構非線性區域規模的大小。式(3)等號左邊括號中表達式的計算結果與結構當前的切向剛度Kt等效，該表達式代表了切向剛度的低秩攝動。當結構局部非線性區域的規模較小時，非線性自由度數目m將遠低于結構的位移自由度數數目n(即m?n)，可利用式(4)對式(3)進行高效求解。

(4)

(5)

當結構規模較大時，雖然式(4)中整體剛度Ke的規模較大，但其為常數矩陣，可僅在分析前合成并進行一次分解，可見，Woodbury公式的使用可以避免傳統非線性分析方法中最為耗時的大規模整體剛度矩陣反復更新和分解運算，取而代之的是僅對小規模Schur補矩陣進行相應運算，實現了計算降維，從而大幅降低非線性分析所需計算消耗。對于地震反應分析問題，需結合運動方程并引入Newmark方法構造動力控制方程，其對應的Woodbury求解公式具有與式(4)完全一致的矩陣表達形式和計算特點。

2 Woodbury+OpenMP并行分析方法

2.1 OpenMP并行機制

OpenMP是適用于共享內存多處理器并行計算的程序設計標準，具有靈活度高、易于編程實現等優點。OpenMP采用Fork-Join的并行模式(Fork即創建或喚醒派生線程，Join即多線程的匯合)，如圖2所示，其基本思路為：程序開始時首先由一個主線程執行，當遇到某個并行任務時，程序將派生出多個線程，此時主線程和派生線程將同步進行程序計算，當該并行任務結束時，派生線程退出工作，主線程繼續執行后續程序計算任務，遇到下一個并行任務后派生線程將被再次喚醒并參與計算，整個計算程序中可以包含多個并行任務。

圖2 OpenMP并行機制

2.2 Woodbury方法的并行實現

結構進入非線性狀態后，需通過迭代完成每個增量計算步的求解。基于Newton-Raphson迭代格式的Woodbury方法非線性分析流程，如圖3所示。由圖3可知，在開始進行迭代求解之前需首先對結構初始彈性剛度矩陣Ke進行LDLT分解，將分解結果進行存儲，以便于在迭代求解過程中隨時調用。在每次迭代時，主要包含3個計算步驟，分別是Woodbury公式計算、單元狀態確定、非線性相關系數矩陣更新，見圖3。本節通過分析上述3個步驟的計算特點，分別設計了基于OpenMP的并行加速方法，各部分具體并行實現方式詳述如下。

圖3 Woodbury方法的迭代求解流程

2.2.1 非線性系數矩陣的并行更新策略

(6)

(7)

對于系數矩陣Kinf，其為對稱矩陣，將式(6)代入式(5)，可得該矩陣的計算表達式

(8)

式中，Aji為與單元j和i有關的子矩陣塊，代表單元j對單元i的影響程度，其計算表達式為

1≤i≤p

(9)

圖4 矩陣Kinf更新過程示意圖

假設在當前迭代計算步中僅單元p為新進入非線性的單元，同時考慮矩陣Kinf的對稱特性，則根據式(8)可知此時僅需計算其第p列數據(包含子矩陣塊A1p，A2p，…，App)并對計算結果進行稀疏近似即可。根據式(9)可知各子矩陣的計算互不干擾，為此，本文采用如下步驟實現這一計算過程的并行化。

步驟1首先計算式(10)。考慮到分析前已完成初始彈性剛度矩陣Ke的分解運算，可直接調用矩陣分解結果通過回代完成式(10)計算，回代過程的計算消耗極小，且可實現基于OpenMP并行加速，本文方法利用Intel MKL中Pardiso求解器實現并行回代。

(10)

式中，Sp為計算Kinf第p列數據時所需中間變量。

步驟2計算對角子矩陣App

(11)

步驟3在程序中創建并行區并添加OpenMP指令對，將子矩陣A1p，A2p，…，A(p-1), p的求解及其稀疏判定過程分發到不同線程上進行并行計算。稀疏判定是通過將每個計算出的子矩陣與式(11)計算結果(即App)進行對比實現的，如圖5所示。圖5中：‖·‖為矩陣范數；τ為預定義的稀疏判定閾值。由圖5可知，若某個計算出的子矩陣Aip(1≤i≤p-1)的范數與相應對角塊App范數之比小于閾值τ，則可將其近似置零。在程序實際執行時，將這些近似置零的子矩陣直接刪除，以降低數據存儲需求和后續Woodbury公式的計算復雜度。

圖5 系數矩陣矩陣Kinf的并行更新過程

當某個迭代步中有多個新進入非線性的單元時，對每個單元的相關子矩陣列均執行上述過程，即可實現系數矩陣Kinf的更新。從論述可知，Woodbury公式中與結構非線性相關的系數矩陣均可劃分為若干相互獨立的子矩陣，各子矩陣在計算時不存在數據交互，因此這些矩陣的更新過程具有較好的可并行特征。上述對于矩陣Kinf的稀疏化處理將導致每次迭代時Woodbury公式計算結果存在一定近似誤差，但該誤差可通過少量增加迭代次數或引入適當的誤差修正策略予以消除，本文采用Yu等提出的基于縮減基的誤差修正方法消除這一近似誤差帶來的不利影響。

本文方法無需在分析前預先設定非線性區域，而是在每次迭代時根據單元狀態確定結果實時確定進入非線性的單元位置，根據論述可知，非線性相關系數矩陣更新過程所需計算量與非線性單元的位置無關，僅與這些單元的數量有關，因此該部分計算效率將隨非線性區域規模的擴大逐漸降低。

2.2.2 Woodbury公式的并行計算

式(4)僅為Woodbury公式的數學表達，為充分發揮其高效性優勢，在實際程序實現時應避免執行耗時的矩陣與矩陣乘法運算。為此，將式(4)從右向左拆解為如表1所示6個計算步，每次迭代時依次執行這6個步驟即可得到Woodbury公式的計算結果。表1中，向量ΔX1～ΔX5均為這一計算過程的中間變量。由表1可知，各計算步僅涉及簡單的矩陣與向量或向量與向量之間的運算，每個計算步的計算量均極小。考慮到表1各計算步中涉及的系數矩陣均具有較高稀疏性，本文采用壓縮存儲格式CSR3(compressed sparse row)對矩陣進行存儲并基于此開展相關矩陣運算，以避免零元素相關運算對計算資源的占用，降低計算復雜度。

表1 Woodbury公式計算步驟

對于表1中的第1個和第5個計算步，可直接調用矩陣Ke的分解結果通過回代計算完成，本文基于Pardiso求解器使用并行回代方法進行這兩個計算步的快速求解；對于表1中的第2個和第4個計算步，其均為稀疏矩陣與向量之間的乘法運算，根據矩陣運算基本原理可知，等號右邊向量ΔX2和ΔX4中各元素的求解過程均互不干擾，因而可基于OpenMP并行方法實現這兩個計算步的加速求解；表1中的第3個計算步要求對以Schur補矩陣為系數矩陣的方程組進行求解，可采用基于OpenMP的并行分解和回代方法進行高效計算，本文方法利用Pardiso求解器實現該計算步的并行求解，由于Schur補矩陣的規模通常較小，且可近似化為具有高度稀疏性的矩陣，該步的計算耗時通常極小；表1中的第6個計算步為兩個向量的加法運算，其計算量極小，且亦可直接實現OpenMP的并行加速。綜上可知，通過將Woodbury公式拆解為6個連續執行的計算步，并基于OpenMP方法對每個計算步的求解進行加速，可實現該公式計算的全過程并行化。考慮到基于矩陣稀疏存儲進行相關運算可大幅降低計算過程的時間和空間復雜度，本文方法可極大提高每次迭代時Woodbury公式的實際程序運行效率。

2.2.3 單元狀態確定過程的并行實現

單元狀態確定的主要作用是在每次迭代時根據計算出的結點位移數據對各單元材料應力、應變、切向模量等參量進行更新，并計算各單元的結點恢復力，以便于更新不平衡力。對于Woodbury方法，還需在單元狀態確定過程中判定各單元是否進入非線性，以便于在后續迭代計算中對Woodbury求解公式中非線性相關矩陣進行更新。由于單元狀態確定過程中各單元之間不存在數據交互，具有典型的并行化特點，因而可直接在該部分程序代碼的開始和結束部位添加OpenMP并行指令對，以實現計算過程的并行加速。

3 算例分析

參考文獻[49]建立如圖6所示的巨型鋼框架結構分析模型，結構共計48層，層高4 m，總高度為192 m，平面尺寸為40 m×40 m，每邊均為5跨，跨度均為8 m，圖6(b)和圖6(c)分別給出了該結構的平立面布置。結構1～24層柱采用Q345鋼，其他層均采用Q235鋼。使用基于隔離非線性法的纖維梁單元建立該結構分析模型，采用考慮隨動硬化的雙線性彈塑性本構模型模擬材料非線性行為，其中屈服后剛度系數設為0.01。為充分驗證本文建立的并行Woodbury方法對于大規模結構非線性問題的高效性，對模型網格進行適當加密，本例中，將各構件每隔1 m劃分一個單元，共計劃分40 608個單元，結構的總位移自由度數超過21萬(共計215 712)。選取El-Centro地震記錄對結構進行非線性地震反應分析，考慮施加平面雙向地震激勵，其中X方向地面運動的峰值加速度調幅至6.0 m/s2，Y方向地面運動的峰值加速度調幅到5.1 m/s2。計算步長設為0.01 s，截取前20 s地震記錄進行分析，共計分析2 000個增量計算步。為驗證本文方法的準確性，同時使用有限元軟件ABAQUS建立該結構的數值分析模型并選用相同地震記錄進行非線性分析，ABAQUS軟件模型的網格密度與本文方法相同，單元類型為B32。

(a) 三維模型

本文Woodbury+OpenMP方法、串行Woodbury法和ABAQUS軟件分析得到的結構頂點沿X和Y方向的位移時程曲線，如圖7所示。由圖7可知，并行和串行Woodbury法的分析結果相同，且兩者與ABAQUS軟件分析結果基本一致。本文方法和ABAQUS軟件算得的結構X方向最大頂點位移分別為0.429 9 m和0.425 3 m，僅相差1.08%，Y方向的結構最大頂點位移分別為0.366 7 m 和0.368 5 m，僅相差0.48%，這表明本文方法計算結果具有與ABAQUS軟件一致的精度水平。本文方法分析時非線性自由度數目的時程曲線，如圖8所示。由于非線性自由度數等于Woodbury求解公式中需實時更新和分解的Schur補矩陣階數，該值越小表明Woodbury公式的求解效率越高，由圖8可知，結構非線性自由度數的峰值為8 379，僅為結構位移自由度數目的3.88%，遠小于傳統非線性分析方法中需實時更新和分解的整體剛度矩陣階數(即位移自由度數)，這說明對于局部非線性問題，使用Woodbury公式能夠在非線性分析時實現大幅度的計算降維，這也是其計算高效性的主要來源。

(a) X方向

圖8 非線性自由度時程曲線

本文方法的所需計算量與結構非線性區域規模有關，對于和本例類似的局部非線性問題，其高效性優勢尤為顯著，雖然隨結構非線性區域規模的提高其計算效率將逐漸降低，但通過結合適當的計算優化策略，如本文采用的Schur補矩陣稀疏優化方法，本文方法的高效性在非線性區域規模較大時依然能夠得到保證。若結構出現全局非線性或大部分單元進入非線性的情況，本文方法的適用性將降低，然而，對于工程結構的抗震分析問題，這一極端情況通常較少出現。

為驗證本文提出的并行加速策略的有效性，本文Woodbury+OpenMP方法、串行Woodbury方法和ABAQUS軟件的總計算耗時，如表2所示。使用的計算平臺為裝配Intel Core i9-10920X處理器的個人PC機，處理器核數為12。可以看出，本文方法僅耗時967 s，為串行Woodbury分析方法計算耗時的18%，這表明對于該算例，本文并行方法的加速比為5.5。此外，ABAQUS軟件的計算耗時約為9.5 h，分別為串行和并行Woodbury方法計算耗時的6.3倍和34.6倍，并行與串行Woodbury方法各部分的計算耗時，如表3所示。可以看出，本文方法在單元狀態確定部分的并行效率最高，這主要是由于該部分計算時各單元完全獨立，不存在串行計算代碼，無需執行額外的進程創建與銷毀工作。由表3可知，系數矩陣更新和Woodbury公式計算的并行效率雖然相對較低，但此兩部分的計算耗時依然能夠實現較大幅度的降低(分別降低了60.4%和49.8%)。由于Woodbury求解公式本身僅需極小的計算消耗，串行分析程序的主要計算耗時集中于單元狀態確定部分和矩陣更新，因而從這兩部分節約出的較高計算時間有助于大幅提升本文方法的整體并行計算效率。

表2 計算耗時對比

表3 Woodbury方法不同計算部分耗時統計

4 結 論

本文利用OpenMP并行編程模式，針對基于Woodbury公式的結構非線性動力反應分析方法，通過對其迭代求解過程中的主要3個主要計算部分(非線性相關系數矩陣更新、Woodbury公式計算和單元狀態確定)進行并行加速，提出了一種用于高效非線性地震反應分析的Woodbury+OpenMP方法。本文具體結論如下：

(2) Woodbury公式在實際應用時可以拆解為依次執行的6個計算步，各計算步僅涉及矩陣與向量或向量與向量之間的運算，由于各計算步均可實現基于OpenMP并行化，本文方法能夠大幅提升該公式的實際程序執行效率。

(3) 本文方法中，Woodbury公式的使用避免了傳統非線性地震反應分析方法所需的整體剛度方程反復更新求解，進一步結合并行計算技術，本文方法在保證計算精度與傳統非線性分析方法相當的前提下可大幅提升結構地震反應分析效率。