Aztec在混凝土細觀數(shù)值模擬中的應用研究

吳建平，趙軍，宋君強，張衛(wèi)民，馬懷發(fā)

1.國防科技大學計算機學院，長沙410073

2.中國水利水電科學研究院流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京100038

1 介紹

混凝土性能研究是混凝土大壩地震破壞機理研究的關(guān)鍵技術(shù)之一，細觀數(shù)值模擬作為一種重要研究手段，將混凝土試件看成由骨料、沙漿及其交界面組成的三相復合材料[1-3]。在模擬試件時，必須對各種材料的分布狀況進行假設(shè)，假設(shè)骨料分布、骨料尺寸等參數(shù)均具有隨機性，即隨機骨料隨機參數(shù)模型已經(jīng)證明很有效[2-3]。該模型以有限元離散為基礎(chǔ)，模擬過程中需要求解對應于全局剛度矩陣的稀疏線性方程組，以算出每一個加載步上的三個位移分量。該線性方程組的求解在總模擬時間中所占比重非常大，已經(jīng)超過90%[4]。

對三維試件數(shù)值模擬中的稀疏線性方程組，采用直接解法時，不僅計算量很大，內(nèi)存需求也很大。而且，為減少計算量和存儲需求，常采用排序技術(shù)來減少矩陣分解過程中的填入元數(shù)量，但即使采用排序技術(shù)，對三維問題，計算量和內(nèi)存需求依然相當龐大，難以接受。與直接解法相反，迭代法只要給定初始解向量和迭代計算規(guī)則，就可以依次不斷計算出新的迭代值，在達到一定的停機準則時就可以停止迭代。迭代過程中，主要的計算量在矩陣向量乘上，這不會改變稀疏矩陣的非零元結(jié)構(gòu)，所以，計算量與系數(shù)矩陣中的非零元個數(shù)成正比，且只需要存儲系數(shù)矩陣與一些相應的向量，從而計算量和存儲需求都相對較少[5]。因此，對三維問題，采用迭代求解是最佳選擇。但迭代法面臨收斂速度慢，甚至不收斂的問題。為此，需要引入預條件技術(shù)，以加速收斂過程。

Aztec是一個迭代求解線性方程組的軟件包[6-7]，該軟件包試圖避開煩瑣的并行編程細節(jié)。在創(chuàng)建分布矩陣、相應的右端向量和解向量，并進行相應的參數(shù)和選項設(shè)置后，就可以調(diào)用Aztec的求解器接口進行求解。Aztec包含多種K rylov子空間迭代法，包括CG、GMRES、BiCGSTAB等，這些方法可以與各種預條件結(jié)合使用，采用的預條件包括多項式預條件，在子區(qū)域內(nèi)采用ILU的區(qū)域分解預條件，在子區(qū)域內(nèi)采用GS迭代的預條件，以及基于最小二乘的預條件等。

雖然Aztec已經(jīng)實現(xiàn)了并行計算，但是該軟件包內(nèi)部并未提供圖分割接口。同時，由于混凝土試件采用的是無結(jié)構(gòu)網(wǎng)格，而已有研究表明，采用M eTis軟件包中基于嵌套二分的圖分割算法最為有效[8]。因此，本文先利用圖分割算法對混凝土試件進行圖分割，之后在該分割下導出相應的分布矩陣和分布向量，相應的矩陣裝配采用基于稀疏數(shù)據(jù)結(jié)構(gòu)的高效并行算法[4]，之后將稀疏矩陣和右端向量導入到Aztec的求解器接口子程序中進行求解。

2 混凝土試件細觀數(shù)值模擬

在混凝土細觀數(shù)值模擬中，將混凝土試樣看成是由骨料、灰漿以及其間界面組成的三相復合材料，這里以隨機骨料隨機參數(shù)模型為基礎(chǔ)[1]。立方體試件用有限元離散，含骨料的區(qū)域用較小的六面體單元進行離散，只含灰漿的區(qū)域用較大的六面體單元進行離散，這兩類離散區(qū)域用四面體單元連接。

所有有限元分為三類：骨料單元、灰漿單元和界面單元，分別對應于材料的三相。對每個有限元，要么有8個形函數(shù)，要么有4個形函數(shù)，這由所對應的單元是六面體單元還是四面體單元來確定。有限單元的每個節(jié)點對應于應力或應變的3個分量。對每個有限元應用損傷模型，并對每個節(jié)點應用簡單的彈性關(guān)系，之后，就可以形成單元剛度矩陣。當每個節(jié)點上的荷載已知之后，節(jié)點位移可以通過求解對應于全局剛度矩陣的稀疏線性方程組來獲得。

2.1 靜載時的數(shù)值模擬

對靜態(tài)加載，荷載一步一步地增加。在加載過程中，如果一個單元所受應力超過其所能承受的極限強度，單元剛度就會退化，即剛度是損傷參數(shù)的函數(shù)。因此，平衡方程是非線性的，采用增量方法進行求解。由于在荷載更大時產(chǎn)生了新的損傷，求解過程需要進行迭代。

用[Ki]、{Ui}、{Pi}分別表示單元剛度矩陣、位移和第i個加載步時的靜荷載，則第i個加載步時的平衡方程可以寫為：

且第i+1個加載步時的平衡方程可以寫為：

因此

這里

且e為單元，e為單元的應變，d是對應于應變的損傷系數(shù)。這樣，第i個加載步的求解過程可描述如下：

2.2 動態(tài)加載時的數(shù)值實驗

在動態(tài)加載時，t時候與t+Dt時候的平衡方程可以分別寫為：

因此

其中

[M]為單元質(zhì)量矩陣，[C]為單元阻尼矩陣，[Kd]為單元動態(tài)剛度矩陣，[Ks]為單元靜態(tài)剛度矩陣，{Pd}為動態(tài)荷載，[Ps]為靜態(tài)荷載。選取δ≥0.5，γ=0.25(0.5+δ)2，并記

則方程（3）可以近似為：

方程（4）是非線性的，可以利用下述迭代進行近似

其中

除了需要求解形如式（5）的線性方程組外，對[M]和[C]的計算還需要計算[K]的某些特征值，這也需要求解相應的線性方程組。

3 Aztec軟件包

Aztec是一個迭代軟件包，用于并行求解稀疏線性方程組

其中，A是n×n稀疏矩陣，b是n維向量，x是要求解的n維向量。Aztec試圖避開煩瑣的并行編程細節(jié)，因此最復雜的環(huán)節(jié)是為特定應用問題進行矩陣分布狀態(tài)的指定。一旦創(chuàng)建了分布矩陣，就可以在任意運行Aztec的并行計算機上進行計算。Aztec提供了兩種稀疏矩陣格式，一種是逐點修正型稀疏行DMSR格式，另一種是逐塊行DVBR格式，分別是MSR格式和VBR格式在分布存儲環(huán)境下的版本。由于現(xiàn)有混凝土細觀數(shù)值模擬程序中采用CSR格式，且CSR格式轉(zhuǎn)換到MSR格式代價很低，因此，文中采用DMSR格式進行分布矩陣指定。

任意長為n的向量通過某種分割方法，將元素分配到特定的處理器上。當計算向量y=Ax中的分量時，一個處理器僅計算分配到其上的y分量，這些分量顯式地存儲在該處理器上，分量對應的全局標號存儲在update索引集中。update又分為兩個子集，即internal和border。集合internal中每個索引對應的分量只需要利用本處理器上的信息進行更新。集合border中含有的索引對應于要在矩陣向量乘中進行更新，但需要從其他處理器獲取x分量值的y分量。在矩陣向量乘中，不在本處理器上，但又需要用來更新border中y分量的x分量指標存儲在索引集external中。這些x分量顯式存儲在其他處理器上，且在進行矩陣向量乘時，需要從這些處理器經(jīng)由通信來獲得。Aztec在對向量進行存儲時，先存儲internal對應的分量，之后存儲border對應的分量，最后存儲external對應的分量。同時，所有從同一個處理器接收到的external分量進行連續(xù)存儲。

Aztec包含許多K rylov子空間迭代法，包括CG、GMRES、BiCGSTAB、CGS和TFQMR。這些K rylov方法可以與各種預條件結(jié)合使用，可用的預條件包括Jacobi迭代，Neumann展開多項式，最小二乘多項式，非重疊區(qū)域分解k步對稱化GS和區(qū)域分解預條件，在區(qū)域分解預條件中，每個子區(qū)域上的局部預條件可以選用LU分解、ILUT、ILU(k)、RILU(k，ω)、BILU(k)和ICC(k)，而對子區(qū)域解到全局解的組合可以采用由擁有者確定最終分量值的標準方式和對分量分別各自進行相加的對稱化方式兩種方式，區(qū)域分解中的重疊度可以采用參數(shù)來進行指定。

4 數(shù)值實驗

針對非結(jié)構(gòu)網(wǎng)格上混凝土細觀數(shù)值模擬所得稀疏線性方程組，就Aztec軟件包中提供的預條件和迭代法組合進行實驗。由于這里的稀疏線性方程組全都對稱正定，因此迭代法全部選用CG法。同時，為利用預條件迭代進行求解，預條件本身必須具有對稱性，因此僅對Aztec中具有對稱性且相對比較有效的預條件進行測試，包括區(qū)域分解預條件，非重疊區(qū)域分解k步對稱化GS即SGS（k）和最小二乘k次多項式即LS（k）。區(qū)域分解預條件中采用對分量各自分別進行相加的對稱化方式，對重疊度0即無重疊的DDM（0）與重疊度1的DDM（1）進行測試，且對局部是否采用RCM排序進行測試，在子區(qū)域局部對具有對稱性的ICC（0）和ICC（1）進行測試。在所有實驗中，初始猜測都選為0向量，且迭代的停機準則選為殘量2范數(shù)下降10-4倍。

所有測試結(jié)果均在一機群系統(tǒng)上得到，該系統(tǒng)采用Intel?Xeon?CPU X 7350@2.93GHz（Cache 4 096 KB），且用Infiniband進行互連，操作系統(tǒng)為Red Hat 4.1.2-44版的Linux，消息傳遞接口采用MPICH2 1.4.0rc2，所用編譯器為Intel Fortran 10.1.018。

實驗針對兩個混凝土試件進行，離散網(wǎng)格形狀參見文獻[3]，試件1是一個三級配混凝土，其大小為1 100mm×300 mm×300mm。離散節(jié)點有44 117個，有限單元有53 200個。兩條支撐軸分別位于x=-0.36m且z=0m，x=0.54m且z=0m處，即處于底部且分別離左右邊界0.1m處。兩條加載軸分別位于x=-0.06m且z=0.3m，x=0.24m且z=0.3m處，即處于試件頂部且分別距左右邊界0.4m。試件2是一濕篩混凝土，其大小為550mm×150 mm×150mm。離散節(jié)點數(shù)為71 013，有限單元數(shù)為78 800。兩支撐軸分別位于x=-0.15m且z=0m，x=0.3m且z=0m處，兩加載軸分別位于x=0m且z=0.15m，x=0.15m且z=0.15m處。

在對兩試件的靜載實驗中，每步增加的荷載為0.25 kN。當在某步出現(xiàn)損傷單元時，在該步上就需要連續(xù)求解多個線性方程組，以校正由剛度矩陣退化引起的位移。對試件1，在第439步出現(xiàn)損傷單元，并在第567步完全失效，總共需要求解696個線性方程組。對試件2，在第59步出現(xiàn)損傷，且在第94步完全失效，總共需要求解178個線性方程組。在對試件2的動載實驗中，時間步長取為0.001 s，且荷載線性地以增量0.8 kN的方式增加，在第38個加載步完全失效，由于實際計算之前，需要計算某些特征值以確定相關(guān)參數(shù)，這也需要求解線性方程組，所以模擬過程總共需要求解221個線性方程組。

對試件1、試件2的靜載實驗，以及對試件2的動載實驗結(jié)果分別列于表1、表2、與表3，其中時間單位為s，“—”表示迭代未能正常結(jié)束。從表中實驗結(jié)果可見，采用RCM排序時，從迭代次數(shù)上看，對LS和SGS預條件，與不采用RCM排序時完全相同。對基于區(qū)域分解的并行ICC預條件，迭代次數(shù)一般稍有減少。從計算時間上看，僅在局部采用ICC（1）且采用重疊區(qū)域分解時才具有改進效果，在其他情況下求解時間一般反而較長，這可能是因為并行預條件的有效性需要全局協(xié)調(diào)一致才能充分體現(xiàn)所致，即使局部預條件很有效，如果區(qū)域分解時策略適當，導致子區(qū)域間對應的舍棄元素相對較大，將對全局計算效果的改進形成很大限制。

此外，從實驗結(jié)果中還可以看到，當處理器個數(shù)增加1倍時，有時出現(xiàn)求解時間減少1倍多，即超線性加速比現(xiàn)象，這應是Cache的影響所致。同時，可以注意到，在很多情況下，在采用相同個數(shù)處理器時，ICC（0）結(jié)合無重疊區(qū)域分解時的求解時間在所有預條件方法中最短，這為以后進行針對混凝土材料細觀數(shù)值模擬時，進行預條件技術(shù)的研究提供了比較基準，即當與Aztec軟件包進行求解性能比較時，只要與其中不采用RCM排序的無重疊區(qū)域分解型并行ICC（0）預條件CG迭代進行比較。

表1 對試件1進行靜載模擬時平均每個線性方程組的迭代次數(shù)和求解時間

表2 對試件2進行靜載模擬時平均每個線性方程組的迭代次數(shù)和求解時間

表3 對試件2進行動載模擬時平均每個線性方程組的迭代次數(shù)和求解時間

5 結(jié)論

本文利用Aztec軟件包中的各種預條件技術(shù)，對目前混凝土細觀數(shù)值模擬中大型對稱正定稀疏線性方程組的求解，進行了實驗比較研究。結(jié)果表明，在基于區(qū)域分解的并行ICC、無重疊對稱化GS迭代、最小二乘等預條件技術(shù)中，采用無重疊區(qū)域分解型并行ICC時求解時間最短，且采用RCM排序一般不僅不能縮短求解時間，反而常導致求解時間更長，因此，在本應用中不采用RCM排序更有利。

[1]唐春安，朱萬成.混凝土損傷與斷裂－數(shù)值模擬[M].北京：科學出版社，2003.

[2]馬懷發(fā)，陳厚群.全級配大壩混凝土動態(tài)損傷破壞機理研究及其細觀力學分析方法[M].北京：中國水利水電出版社，2008.

[3]馬懷發(fā).全級配大壩混凝土動態(tài)性能細觀力學分析研究[D].北京：中國水利水電科學研究院，2005.

[4]吳建平，王正華，朱星明，等.混凝土細觀力學分析程序中的快速算法與并行算法設(shè)計[J].計算力學學報，2008，25（3）：352-358.

[5]吳建平，王正華，李曉梅.稀疏線性方程組的高效求解與并行計算[M].長沙：湖南科技出版社，2004.

[6]Wu Jianping，Zhao Jun，Song Junqiang，et al.Impact of two factors on several domain decomposition based parallel incomplete factorizations for the meso-scale simulation of concrete[C]//Proceedings of the 3rd International Conference on Information and Computing Science（ICIC2010），Wuxi，China，2010.

[7]Shadid J N，Tum inaro R S.Aztec-a parallel preconditioned K rylov solver library：algorithm description version 1.0，Technical Report Sand95[R].Sandia National Laboratories，Albuquerque NM，87185，1995.

[8]Tum inaro R S，Heroux M，Hutchinson S A，et al.Official Aztec user’s guide，Version 2.1，Technical Report Sand99-8801J[R].Sandia National Laboratories，Albuquerque NM，87185，1999.