徐鳳洲，張健飛

(河海大學 力學與材料學院，江蘇 南京 211100)

0 引 言

材料和結(jié)構(gòu)破壞的計算機數(shù)值模擬一直是傳統(tǒng)數(shù)值方法面臨的一個難題，其局限性源自連續(xù)介質(zhì)力學理論中的連續(xù)性假設和實際問題不連續(xù)之間的矛盾，在分析裂紋問題時，裂紋的起裂和擴展都必須引入外部的準則進行判斷，無疑增加了分析問題的復雜性。近場動力學(peridynamics，PD)作為一種新興的基于非局部作用的理論[1]和方法，在分析不連續(xù)問題時，有其獨特的優(yōu)勢。近場動力學的運動(平衡)方程是采用空間積分形式，并基于牛頓第二定律，因此，對連續(xù)體和非連續(xù)體均適用。因其本構(gòu)模型中已經(jīng)包含了對損傷和斷裂的描述，所以不需要附加開裂判據(jù)和裂紋路徑分析，裂紋自然而然地萌生并擴展[2]。然而，近場動力學在分析裂紋破壞等不連續(xù)問題的過程中，需要將模型離散為一系列包含物質(zhì)信息的物質(zhì)點，進一步需要搜索每個物質(zhì)點鄰域范圍內(nèi)的，與其相互作用的其他物質(zhì)點，同時，在求解計算過程中，為了保證數(shù)值準確性和穩(wěn)定性還需要進行大量的迭代求解[3]。因此，近場動力學在實際數(shù)值模擬中面臨著計算復雜度高、計算量大等問題[3-4]，這無疑會有礙于近場動力學的進一步發(fā)展和實際應用。

并行計算是目前擴大計算規(guī)模和加快計算速度的重要途經(jīng)，開展PD數(shù)值模擬的并行計算研究可以有效擴大計算規(guī)模和提高計算速度。目前，PD數(shù)值模擬的并行計算研究采用的方法主要有MPI并行[5]、基于GPU和 CUDA的多線程并行[3]和基于GPU和OpenCL的多線程并行[6]等，但是這些方法的實現(xiàn)過程較為復雜，且對計算設備要求較高。OpenMP(open multi-processing)是一套用于共享內(nèi)存并行系統(tǒng)多線程程序設計的指導性注釋(compiler directive)，是為共享內(nèi)存的多處理器系統(tǒng)設計的并行編程方法。OpenMP特別適合用在多核處理器計算機上運行的并行程序設計，可以快速實現(xiàn)PD模擬的并行化。目前，基于OpenMP的PD并行計算研究較少，缺乏具體實現(xiàn)過程描述。本文針對PD計算特點，對PD程序結(jié)構(gòu)做歸納和總結(jié)，并基于OpenMP對PD數(shù)值模擬的并行實現(xiàn)方法進行研究，對程序的主要模塊做并行化，編制并行程序，通過經(jīng)典算例測試程序的性能，以期能提高計算效率。

1 基于鍵型的PD串行算法

近場動力學理論目前主要分為鍵型PD理論[7]和態(tài)型PD理論[8]?？紤]到鍵型理論較早提出[7]，且研究比較深入和成熟[2]，本文選擇鍵型PD理論，并采用鍵型PMB模型。PD理論的計算框架[9]主要包括模型離散、數(shù)據(jù)輸入和存儲、初始和邊界條件施加、時間積分、物質(zhì)點內(nèi)力求解、破壞模型、數(shù)據(jù)輸出等。另外，模型離散后需要根據(jù)近場尺寸范圍搜尋并構(gòu)建鄰接節(jié)點域[9]。根據(jù)已有的PD計算框架，本文首先采用C++程序語言編制鍵型PD串行程序，并采用準靜態(tài)模擬和動態(tài)裂紋擴展模擬2種模型。前者適用于不考慮損傷情況的一般彈性變形問題，后者則適用于快速加載條件下的裂紋擴展問題。以動態(tài)裂紋擴展模擬為例，PD程序的大體結(jié)構(gòu)如圖1所示。

圖1 PD程序結(jié)構(gòu)圖

1.1 準靜態(tài)模擬

對于準靜態(tài)模擬，本文采用自適應動態(tài)松弛法求解問題，不考慮物質(zhì)點的損傷情況，適用于模擬一般彈性變形過程，例如板的軸向受拉問題等。對于準靜態(tài)問題，PD的運動(平衡)方程為

(1)

在準靜態(tài)模擬中，程序按執(zhí)行順序大體可以分為以下幾個模塊：

(1)讀取模型信息，將模型離散化——Data input模塊。

(2)構(gòu)建每個物質(zhì)點的鄰接節(jié)點域——Neighbor模塊。

(3)計算每個物質(zhì)點的表面修正系數(shù)，改善PD模型自身的邊界效應問題——Compute kave模塊。

(4)參數(shù)初始化——Initial模塊。

(5)時間積分，更新物質(zhì)點位置。其中又包含邊界條件施加、計算物質(zhì)點內(nèi)力密度和自適應阻尼項、檢查收斂性、更新物質(zhì)點位移值，結(jié)果輸出等子模塊——Time integration模塊。

1.2 動態(tài)裂紋擴展模擬

在模擬裂紋擴展問題中，需要考慮物質(zhì)點的損傷情況，主要有2種類型：一種是準靜態(tài)裂紋擴展；一種是動態(tài)裂紋擴展。本文采用動態(tài)裂紋擴展模擬，這是因為動態(tài)裂紋擴展應用較為廣泛，與準靜態(tài)模擬不同，它是采用顯式求解方法，并且裂紋分叉結(jié)果更加明顯、更加直觀、可視化更好。對于動態(tài)問題，PD的運動(平衡)方程為

(2)

式中各參數(shù)的意義與準靜態(tài)問題相同。在動態(tài)裂紋擴展模擬中，程序的大體結(jié)構(gòu)與上述準靜態(tài)模擬相似，但由于求解方法不同，在時間積分模塊(Time integration)中，去掉了計算自適應阻尼項、檢查收斂性等子模塊，并且在計算物質(zhì)點內(nèi)力密度子模塊中，增加了表征物質(zhì)點損傷的量，需要分別對鍵破壞以及物質(zhì)點損傷作出判斷和計算。

2 OpenMP并行編程模型

OpenMP是一種面向共享內(nèi)存以及分布式共享內(nèi)存的、多處理器多線程并行語言。其編程模型是以線程為基礎，通過編譯制導語句在原來串行程序的基礎上顯示進行并行化處理[10]。對原有串行程序的改動較小，具有操作簡便、開發(fā)效率高等優(yōu)點。OpenMP采用的執(zhí)行模式是串行→并行→串行→并行→串行……，如圖2所示。這種執(zhí)行模式的核心在于并行區(qū)域中線程的派生/縮并(Fork/Join)[11-12]。

圖2 OpenMP并行程序的運行過程

3 基于OpenMP的PD并行算法

為了對已編制的PD串行程序進行更好的并行化處理，本文首先測試了PD串行程序各模塊的計算時間，再對計算耗時長的模塊進行并行化。測試結(jié)果顯示,計算耗時長的模塊主要集中在搜索每個物質(zhì)點近場范圍內(nèi)的鄰域節(jié)點(Neighbor模塊)、計算每個物質(zhì)點的表面修正系數(shù)(Compute kave模塊)，以及時間積分(Time integration模塊)這3個模塊，并進一步考查了單個時間步內(nèi)，時間積分模塊(Time integration)中各子模塊的計算時間，結(jié)果發(fā)現(xiàn)，計算物質(zhì)點內(nèi)力密度子模塊的計算耗時占了這一模塊(Time integration)計算耗時的大部分(80%左右，計算設備的性能不同，結(jié)果會存在差異)，所以在對程序進行并行處理時，需要著重對這一部分進行并行優(yōu)化。分析了程序的熱點循環(huán)和熱點函數(shù)后，在對應的模塊內(nèi)，主要針對無循環(huán)依賴的循環(huán)體，增添OpenMP編譯制導語句，對程序進行并行化。

3.1 對Neighbor模塊的并行化

程序中Neighbor模塊的主要作用是搜索每個物質(zhì)點近場范圍內(nèi)的，與其相互作用的其他物質(zhì)點，構(gòu)建鄰接節(jié)點域。在程序的這一模塊中，將構(gòu)建3個數(shù)組分別存放所有物質(zhì)點的鄰域節(jié)點號(按物質(zhì)點排列順序存放)(數(shù)組1)(Neighbor)、所有物質(zhì)點的鄰域節(jié)點個數(shù)(數(shù)組2)(numNeighbor)，以及所有物質(zhì)點的起始鄰域節(jié)點在數(shù)組1中的下標位置(數(shù)組3)(Nodefam)。由于數(shù)組1是存放所有物質(zhì)點的鄰域節(jié)點號，因此,必須按照物質(zhì)點的排列順序，依次將每個物質(zhì)點的鄰域節(jié)點號順序存放到數(shù)組1中。數(shù)組3中存放的是所有物質(zhì)點的起始鄰域節(jié)點在數(shù)組1中的下標位置，由于數(shù)組1是按照物質(zhì)點排列順序存放的，因而數(shù)組3也需要按照順序計算和存放數(shù)據(jù)。

整個Neighbor模塊將由兩個大的循環(huán)體構(gòu)成。第一個循環(huán)體是二重嵌套循環(huán)，主要完成以下過程：搜索每個物質(zhì)點的鄰域內(nèi)節(jié)點，構(gòu)建鄰接節(jié)點域；確定每個物質(zhì)點的鄰域內(nèi)節(jié)點個數(shù)并存放到數(shù)組2中。由于第一個循環(huán)體沒有循環(huán)依賴，每個物質(zhì)點可以獨立運行，因而可以進行并行化，并行化主要過程如圖3所示。第二個循環(huán)體構(gòu)建了上述的數(shù)組1和數(shù)組3，由于數(shù)組1和數(shù)組3中數(shù)據(jù)的存放需要按照物質(zhì)點排列順序，即循環(huán)指標的先后次序，排隊執(zhí)行，因而無法進行并行化，只能依照原有的串行方式執(zhí)行。好在第二個循環(huán)體的運行時間相比第一個循環(huán)體要小得多，一般只有零點幾秒，所以并沒有影響到整個Neighbor模塊的并行優(yōu)化效果，這一模塊的部分并行代碼段如圖4所示。private()括號內(nèi)為聲明的私有變量列表，需填上循環(huán)體內(nèi)的私有變量(特別是進行寫操作的變量)，確保每個物質(zhì)點在程序執(zhí)行過程中互不干擾，避免出現(xiàn)數(shù)據(jù)競爭現(xiàn)象。并行語句for指令后不帶schedule子句，系統(tǒng)會默認采用靜態(tài)調(diào)度方式(調(diào)度過程見圖5)，即如果將循環(huán)迭代的總次數(shù)設為n，并行區(qū)域內(nèi)線程總數(shù)為p，那么分配給每個線程約n/p次連續(xù)的迭代[13]。因此不同處理器在對共享數(shù)組(例如數(shù)組2)進行寫操作時，操作的是數(shù)組中相隔比較遠的元素，從而避免了偽共享[13]的問題。

圖3 Neighbor模塊并行化示意圖(以3個線程為例)

圖4 并行代碼段1(構(gòu)建物質(zhì)點鄰接節(jié)點域)

3.2 對Compute kave 模塊的并行化

Compute kave模塊的主要作用是改善PD模型自身的邊界效應，對物質(zhì)點的微模量系數(shù)值作出相應的修正。為了保證修正系數(shù)的準確性，將這一部分迭代了20次[14]。每一次迭代，每個物質(zhì)點需要對鄰域內(nèi)節(jié)點做一次循環(huán)，所以是一個二重嵌套循環(huán)。OpenMP可以對嵌套循環(huán)體內(nèi)的任意一個循環(huán)進行并行化，為了避免過多的并行開銷導致并行效率降低，本文選擇對外層循環(huán)做并行優(yōu)化，而內(nèi)層循環(huán)依然保持串行執(zhí)行方式。在并行語句for指令后增添schedule(guided)子句，系統(tǒng)將會采用指導性調(diào)度方式進行調(diào)度(調(diào)度過程見圖6)。指導性調(diào)度方式是一種動態(tài)方式，與靜態(tài)調(diào)度方式(schedule(static))不同，分配給每個線程的迭代塊不是均等的，而是在開始時迭代塊會比較大，后來逐漸減小，從而有效地減少了調(diào)度開銷[13]。這樣處理是考慮到每個物質(zhì)點的鄰域內(nèi)節(jié)點個數(shù)存在差異，特別是邊界處的物質(zhì)點，其鄰域內(nèi)節(jié)點個數(shù)會大大減少，如果仍然采用靜態(tài)調(diào)度方式，則難以保證線程在程序執(zhí)行過程中的負載平衡，而動態(tài)調(diào)度方式(schedule(dynamic))又存在額外的開銷，不利于提高并行效率，所以最終考慮采用指導性調(diào)度方式(schedule(guided))。如果需要連續(xù)對多個for循環(huán)體進行并行化，那么在并行語句#pragma omp parallel下方需要增添一對大括號將整個并行區(qū)域括起來。

3.3 對Time integration模塊的并行化

Time integration模塊一般是整個程序中計算時間最長的部分(具體運行時間與設置的迭代步數(shù)有關)。這一模塊的主體是一個大的迭代循環(huán)，可能包含成千上萬次的迭代。迭代過程中，前后兩次迭代是存在相互依賴關系的，并且迭代循環(huán)的次數(shù)在某些情況下也是無法提前確定的，所以對迭代循環(huán)整體無法運用OpenMP進行并行優(yōu)化，但是在每一次迭代時間步下，程序的執(zhí)行過程是可以并行化的。本文主要對Time integration模塊中的計算物質(zhì)點內(nèi)力密度，更新物質(zhì)點位移和速度等子模塊進行并行優(yōu)化，在程序的每一次迭代中這些子模塊都會相應地被執(zhí)行一次。計算物質(zhì)點內(nèi)力密度子模塊的部分并行代碼段如圖7所示。對涉及物質(zhì)點鄰域節(jié)點循環(huán)的for循環(huán)體，本文都是采用指導性調(diào)度方式，而對于一般的物質(zhì)點循環(huán)，例如更新物質(zhì)點位移、速度等都是采用靜態(tài)調(diào)度方式。

圖7 并行代段碼2(計算物質(zhì)點內(nèi)力密度)

在整個程序的并行化和調(diào)試過程中會用到一些OpenMP運行環(huán)境操作函數(shù)，例如omp_set_num_threads(n)，此函數(shù)用來確定并行區(qū)域內(nèi)運行的子線程數(shù)量，只能在并行區(qū)域之外調(diào)用[13]。實際操作中，一般將此函數(shù)放在編譯制導語句#pragma omp parallel的前面，通過設定n值確定并行區(qū)域運行時子線程的數(shù)量。

4 計算結(jié)果與分析

為了更好地考查程序的并行性能，以及計算規(guī)模與并行加速比之間的關系，選擇以下3個算例，物質(zhì)點總數(shù)分別為10 600，103 000，253 000。通過對比串行程序與并行程序的運行計算結(jié)果，驗證并行程序的正確性，計算并分析3個算例對應的并行加速比情況。下面將詳細闡述這3個算例的計算結(jié)果。

4.1 準靜態(tài)模擬計算結(jié)果

第一個算例為準靜態(tài)鍵型PD模型算例，計算對象為二維各向同性軸向受拉板(圖8)，無預制裂紋，板長和板寬均為50 mm，板厚設為單位1，密度ρ=8 000 kg/m3，彈性模量E=192 GPa，泊松比v=1/3。離散化參數(shù)：x方向和y方向的物質(zhì)點離散間距均為0.5 mm，近場范圍為3倍的物質(zhì)點間距。邊界條件：在板的上下兩端虛擬邊界層處施加大小分別為2 000 mm/s和 -2 000 mm/s的速度邊界條件，最大拉伸位移值d=0.2 mm，虛擬邊界層厚度為3倍的物質(zhì)點間距?？偟奈镔|(zhì)點數(shù)為10 600，最大迭代步數(shù)設為5 000，時間步長dt=1.0×10-7s/step。為了模擬板在上下兩端受拉情況下的彈性變形過程，不考慮物質(zhì)點的損傷情況，將臨界伸長率s0設為1。

圖8 二維軸向受拉板

為了驗證程序計算結(jié)果的正確性，分別取y坐標為25.25 mm處的一排點x方向和y方向的位移值，和理論解作對比，對比結(jié)果如圖9所示。

圖9 部分物質(zhì)點位移值PD解和理論解的對比

從圖9可以看出，除了靠近邊界處的物質(zhì)點位移值有微小的差異(誤差<5%)，這是由PD數(shù)值模型自身所引起的，其余物質(zhì)點位移值的PD解和理論解基本保持一致。表1給出了圖8中標注的部分物質(zhì)點在模型達到最終平衡狀態(tài)時，x方向和y方向的位移值dispx，dispy，并將串行程序和并行程序的計算結(jié)果作對比，可以發(fā)現(xiàn)兩者基本保持一致，從而驗證了并行程序的正確性。

表1 串行程序和并行程序計算結(jié)果的對比

4.2 動態(tài)裂紋擴展模擬計算結(jié)果

第二個算例為二維各向同性矩形板裂紋分叉模擬算例。板長和板寬分別為100,40 mm，在板左側(cè)中心處有一條長度lc=50 mm的預制裂紋。板的上下兩側(cè)有均勻拉伸荷載σ=12 MPa，如圖10所示。材料參數(shù)：彈性模量E=65 GPa，泊松比v=1/3，密度ρ=2 235 kg/m3。離散化參數(shù)：x方向和y方向的物質(zhì)點離散間距均為0.2 mm，近場范圍為4倍的物質(zhì)點間距，總物質(zhì)點數(shù)為103 000。時間步長dt=5.0×10-8s/step，迭代步數(shù)設為1 000，臨界伸長率s0=0.001 6。

圖10 軸向拉伸下的單邊預制裂紋二維矩形板

圖11給出了模型在迭代步數(shù)為500和1 000時的裂紋擴展情況?？梢钥闯觯鸭y從預制裂紋的右端開始逐漸延伸，到一定階段開始分叉。由于預制裂紋位置的對稱性，裂紋擴展路徑也呈現(xiàn)出完美的對稱性(與已有文獻[6]圖8中的裂紋分叉結(jié)果基本一致)。

圖11 算例二的計算結(jié)果

第三個算例為板中心裂紋動態(tài)擴展模擬算例，是在準靜態(tài)模擬算例模型的基礎上，在板的中心位置增添了1條預制中心裂紋，如圖12所示，裂紋長度2a=10 mm。離散化參數(shù)：物質(zhì)點離散間距為0.1 mm，近場范圍為3.015倍的物質(zhì)點間距，總物質(zhì)點數(shù)為253 000。邊界條件：在板的上下兩端虛擬邊界層處施加大小分別為50 000，-50 000 mm/s的速度邊界條件。時間步長dt=1.3367×10-8s/step，臨界伸長率s0=0.0447 2，迭代步數(shù)設為1 000，其他參數(shù)保持不變。

圖12 軸向拉伸下含中心裂紋的二維板

表2給出了算例三在線程數(shù)設為1的情況下程序各模塊的運行時間情況。

表2 算例三在線程數(shù)為1時程序各模塊的運行時間

Insert center crack模塊的作用是在板中心位置處添加預制裂紋，其余各模塊的含義在前文已經(jīng)介紹過，在此不再復述。由表2可以看出，在程序的整體運行過程中，相比其他模塊，計算時間比較長的3個模塊分別是Neighbor、Time integration和Compute kave模塊。

圖13給出了模型的初始狀態(tài)，以及迭代步數(shù)為650和1 000時的裂紋擴展情況。從圖13可以看出，裂紋從預制裂紋的兩端開始逐漸擴展、分叉，形狀猶如魚的骨頭(與文獻[15]圖11(b)、文獻[16]圖9.7中的裂紋分叉結(jié)果基本保持一致)，直到板完全被撕裂，整個動態(tài)裂紋擴展過程結(jié)束。

4.3 并行計算結(jié)果與分析

本文所采用的實驗平臺為Intel(R)Xeon(R)Core44 CPU E5-2699 v4，主頻為2.20 GHz；操作系統(tǒng)為Windows10；開發(fā)工具為VS2012(支持OpenMP 2.0)。

上述3個算例程序整體運行時間與線程數(shù)之間的關系如表3所示，程序主要模塊在不同線程數(shù)下的并行加速比情況如圖14所示。

由表3、圖14分析可知：隨著線程數(shù)增加，程序運行時間呈下降趨勢，且計算規(guī)模越大，粒子數(shù)越多，下降幅度越明顯，得到的加速比越大，并行效率越高。在計算規(guī)模很小，物質(zhì)點數(shù)不多的情況下，隨著線程數(shù)增多，并行加速比并沒有得到顯著提高，而是趨于平緩，如圖14(d)算例一所示。

圖13 算例三的計算結(jié)果

表3 不同線程數(shù)下3個算例程序整體運行時間

這是因為運用OpenMP對程序進行并行化的過程中，線程組的派生和縮并，子線程的交互運行都需要花費額外的時間，這些都屬于OpenMP自身的并行開銷，而且隨著線程數(shù)的增加，并行開銷所占比例會不斷增大，因此，當計算規(guī)模維持不變，程序運算負載不夠大時，這些開銷會影響到并行程序的加速執(zhí)行，從而并行效率會降低。所以在問題規(guī)模固定的情況下，增加線程并不能持續(xù)提高計算效率，只有選擇合適的線程數(shù)才能獲得最快的計算速度[17]。因此，在實際并行化操作中，應根據(jù)具體的計算規(guī)模選擇合適的線程數(shù)，避免因線程數(shù)量過多而速度更慢所導致的計算資源浪費和計算時間延長，或者計算規(guī)模應該隨著計算線程數(shù)的增加而相應擴大。

從圖14還可以看出，Neighbor模塊的并行加速比情況比較理想，而其余幾個模塊的并行加速比沒有顯著提高。這是因為OpenMP主要針對無循環(huán)依賴的循環(huán)體進行并行化，如果程序中包含沒有循環(huán)依賴的循環(huán)體，使用OpenMP能夠大幅度地提高程序的并行性能。Neighbor模塊中的循環(huán)體無循環(huán)依賴，而其余幾個模塊不僅有循環(huán)依賴的迭代循環(huán)，而且還存在子模塊的調(diào)用，會增加額外的調(diào)用開銷，因此，Neighbor模塊的并行加速比比其他模塊更為理想。

圖14 3個算例程序中主要模塊在不同線程數(shù)下的并行加速比曲線

5 結(jié) 語

在多核處理器平臺上，基于OpenMP的近場動力學模擬相比串行平臺，效率可以得到較好提升，隨著線程數(shù)量的增加，加速比不斷增大，計算時間大幅減少。這有利于充分利用計算機CPU資源，避免大量的CPU資源被閑置和浪費，也有利于近場動力學理論的進一步發(fā)展和應用。本文針對平面問題的鍵型PD模擬開展研究，今后需要進一步研究空間問題和態(tài)型PD模擬，在大型并行計算機上對并行算法和程序的可擴展性做進一步的測試與優(yōu)化。