王星馳, 黃金晶, 蔣 濤

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

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

1 問題描述

考慮三維對流擴散方程

初始條件為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ù);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個粒子對物理量f及其一階導(dǎo)數(shù)f在空間位置r=(x,y,z)處的粒子近似式[1]為(f)i,α=∑jVj(fj-fi)i,αWij, 其中Vj表示第j個粒子對應(yīng)的體積;i,αWij=?Wij/?xi,α, 式中Wij=W(|ri-rj|,h)為核函數(shù),xi,α表示第i個粒子在位置r處的第α個分量,h為光滑長度, 它決定了W的支持域和影響域．

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

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

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

3 主要結(jié)果

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

3.1 有解析解的三維算例

表1 在中心點(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時, 不同時刻沿對角線x=y的變化曲線Fig.3 Change curve along the diagonal line x=y in different time for ν=0.1

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

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

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

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

3.2 無解析解的三維算例

考慮求解域為Ω=[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

其中源項為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, 時間步長dt=10-3,采用36個CPU．

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

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