王星馳, 黃金晶, 蔣 濤

(揚(yáng)州大學(xué)數(shù)學(xué)科學(xué)學(xué)院, 江蘇 揚(yáng)州 225002)

對流擴(kuò)散方程是偏微分方程的一個(gè)重要分支,在流體力學(xué)[1]、生物醫(yī)學(xué)[2]等領(lǐng)域有廣泛應(yīng)用,故對其數(shù)值計(jì)算方法的研究具有重要的理論和實(shí)際意義．基于網(wǎng)格類的方法已被廣泛應(yīng)用于此類方程的數(shù)值求解[3-4], 但其依賴于網(wǎng)格剖分的特性, 使其對高維問題的數(shù)值模擬比較困難,特別對非規(guī)則區(qū)域問題或區(qū)域上存在不均勻布點(diǎn)的問題[5]．SPH[6-7](smoothed particle hydrodynamics) 作為一種純無網(wǎng)格方法, 具有: 1) 易處理區(qū)域上不均勻布點(diǎn)情況或高維復(fù)雜區(qū)域的問題; 2) 算法易實(shí)現(xiàn); 3) 具有分子動(dòng)力學(xué)方法的相鄰粒子并行搜索特性, 故該算法針對三維問題模擬時(shí), 易采用基于MPI (multi-point interface)的并行技術(shù)來提高計(jì)算效率．本文基于核梯度, 擬通過改進(jìn)傳統(tǒng)SPH方法, 并結(jié)合MPI并行技術(shù), 對三維對流擴(kuò)散方程進(jìn)行高效數(shù)值計(jì)算．

1 問題描述

考慮三維對流擴(kuò)散方程

初始條件為u(x,y,z,0)=φ(x,y,z), (x,y,z)∈Ω; 邊界條件為u(x,y,z,t)=φ(x,y,z,t), (x,y,z,t)∈?Ω×(0,T], 其中t為時(shí)間,ν為黏度系數(shù);bx,by,bz為對流系數(shù), 假設(shè)bi<0(i=x,y,z);f,φ,φ為充分光滑的已知函數(shù);Ω為長方體區(qū)域．

2 本文方法

2.1 CSPH-3D方法

為提高傳統(tǒng)SPH方法的精度, 基于Taylor級數(shù)展開給出了一階修正SPH方法[8]．在傳統(tǒng)SPH方法的數(shù)值模擬中, 第i個(gè)粒子對物理量f及其一階導(dǎo)數(shù)f在空間位置r=(x,y,z)處的粒子近似式[1]為(f)i,α=∑jVj(fj-fi)i,αWij, 其中Vj表示第j個(gè)粒子對應(yīng)的體積;i,αWij=?Wij/?xi,α, 式中Wij=W(|ri-rj|,h)為核函數(shù),xi,α表示第i個(gè)粒子在位置r處的第α個(gè)分量,h為光滑長度, 它決定了W的支持域和影響域．

2.2 CSPH-3D方法的粒子搜索并行計(jì)算

圖1 CSPH-3D方法的流程圖Fig.1 Flow chart of CSPH-3D method

在傳統(tǒng)SPH方法模擬中,相鄰粒子搜索在高維問題下會占用較大計(jì)算內(nèi)存和較長CPU計(jì)算時(shí)間,因此本文根據(jù)模擬方程中粒子位置不動(dòng)的特點(diǎn),程序模擬物理量更新前對每個(gè)粒子的相鄰粒子進(jìn)行了標(biāo)定,隨后結(jié)合MPI并行技術(shù)[10],先對所有粒子進(jìn)行編號,再根據(jù)CPU數(shù)量將粒子平均分配,最后同時(shí)進(jìn)行計(jì)算,以此提高CSPH-3D方法模擬三維對流擴(kuò)散方程的計(jì)算效率．該方法詳細(xì)流程圖如圖1所示．為了提高計(jì)算效率,分別在計(jì)算質(zhì)量、密度、一階導(dǎo)數(shù)及最后計(jì)算函數(shù)值的更新中進(jìn)行了并行優(yōu)化．

3 主要結(jié)果

本文并行計(jì)算采用Red Hat Enterprise Linux5.8x86_64操作系統(tǒng), MPI類型為Intel MPI Toolkits 4.0.3, 共有17個(gè)IBM BladeCenter HS22計(jì)算節(jié)點(diǎn), 每個(gè)節(jié)點(diǎn)包括12個(gè)主頻為2.67 GHz的Xeon 6C X5650 CPU．

3.1 有解析解的三維算例

表1 在中心點(diǎn)(0.5,0.5,0.5)處數(shù)值解與精確解的絕對誤差Tab.1 The absolute error between the numerical solution and the exact solution at the center point (0.5,0.5,0.5)

圖2 不同黏度系數(shù)在沿x=y的變化曲線Fig.2 Change curve along the diagonal line x=y in different viscosity coefficients

圖3 黏度系數(shù)ν=0.1時(shí), 不同時(shí)刻沿對角線x=y的變化曲線Fig.3 Change curve along the diagonal line x=y in different time for ν=0.1

表2給出了在813個(gè)粒子數(shù)下不同CPU數(shù)的消耗時(shí)間．從表2可以看出, 三維的串行程序占用了較大的內(nèi)存和較長的計(jì)算時(shí)間,尤其在第一步粒子搜索時(shí); 而隨著CPU的增加,計(jì)算耗時(shí)迅速減少,計(jì)算效率得到提高．定義: 相對加速比τ=單節(jié)點(diǎn)上并行算法的消耗時(shí)間/N個(gè)節(jié)點(diǎn)上并行算法的消耗時(shí)間,本文每個(gè)節(jié)點(diǎn)上有2個(gè)CPU．結(jié)果表明, 相對加速比低于CPU個(gè)數(shù)增加的比值,這是因?yàn)殡S著CPU個(gè)數(shù)增加,從不同CPU上得到結(jié)果的通訊時(shí)間相對增加．

表3給出了不同粒子數(shù)、不同CPU數(shù)下的平均消耗時(shí)間．從表3中數(shù)據(jù)知, 當(dāng)CPU數(shù)不變時(shí), 隨粒子數(shù)的增加,平均消耗時(shí)間也相對增加; 因此通過增加CPU個(gè)數(shù)進(jìn)行并行計(jì)算可以減少計(jì)算機(jī)的模擬時(shí)間, 提高計(jì)算效率．

表2 粒子數(shù)為813時(shí)不同CPU數(shù)下不同計(jì)算步數(shù)的消耗時(shí)間Tab.2 Consumption time of different computational steps at different CPU number when particle number is 813 s

表3 不同粒子數(shù)、不同CPU數(shù)下的平均消耗時(shí)間Tab.3 Average consumption time at different particle number and CPU number s

3.2 無解析解的三維算例

考慮求解域?yàn)棣?[0,1]×[0,1]×[0,1]的無解析解算例

圖4 隨粒子數(shù)增加沿z軸(x=y=0.5)的變化曲線Fig.4 Change curve along z axis(x=y=0.5) with the increase of particle number

其中源項(xiàng)為f(x,y,z,t)=π[2cos(2πt)sin(πx)+cos(πx)sin(2πt)]sin(πy)sin(πz) +π[cos(πy)·sin(πz)+cos(πz)sin(πy)]sin(2πt)sin(πx)+3vπ2sin(2πt)sin(πx)sin(πy)sin(πz), 黏度系數(shù)ν=10-3, 時(shí)間步長dt=10-3,采用36個(gè)CPU．

圖4分別選取粒子數(shù)為513,713,1013,在t=0.25時(shí)沿z軸的數(shù)值收斂圖, 可以看到除兩端梯度有較小誤差外均擬合較好．圖5給出了粒子數(shù)為513時(shí)z=0.5的截面等值線．從圖5可以看出,不同時(shí)刻截面的變化趨于一致．圖4和圖5結(jié)果表明,本文給出的方法隨時(shí)間增長, 數(shù)值結(jié)果均有較好的收斂性和穩(wěn)定性．

圖5 513粒子數(shù)下不同時(shí)間的截面等值線圖Fig.5 Cross-sectional isolines at different times for 513 particles