蔣濤 黃金晶 陸林廣 任金蓮

(揚州大學數(shù)學科學學院,水利與能源動力工程學院,揚州 225002)

為提高傳統(tǒng)光滑粒子動力學(SPH)方法求解高維非線性薛定諤(nonlinear Schr?dinger/Gross-Pitaevskii equation,NLS/GP)方程的數(shù)值精度和計算效率,本文首先基于高階時間分裂思想將非線性薛定諤方程分解成線性導數(shù)項和非線性項,其次拓展一階對稱SPH方法對復數(shù)域上線性導數(shù)部分進行顯式求解,最后引入MPI并行技術,結合邊界施加虛粒子方法給出一種能夠準確、高效地求解高維NLS/GP方程的高階分裂修正并行SPH方法.數(shù)值模擬中,首先對帶有周期性和Dirichlet邊界條件的NLS方程進行求解,并與解析解做對比,準確地得到了周期邊界下孤立波的奇異性,且對提出方法的數(shù)值精度、收斂速度和計算效率進行了分析;隨后,運用給出的高階分裂粒子方法對復雜二維和三維NLS/GP問題進行了數(shù)值預測,并與其他數(shù)值結果進行比較,準確地展現(xiàn)了非線性孤立波傳播中的奇異現(xiàn)象和玻色-愛因斯坦凝聚態(tài)中帶外旋轉項的量子渦旋變化過程.

1 引 言

眾所周知,非線性薛定諤方程(nonlinear Schr?dinger/Gross-Pitaevskii equation,NLS/GP)常被用來描述一些非線性物理特性,例如：晶體中熱脈沖、非線性光學、等離子和量子動力學現(xiàn)象等[1-5].NLS/GP方程的研究涉及復雜非線性項和外部旋轉導數(shù)項,導致孤立波值隨時間變化過程在許多情況下難以用理論手段[6-8]精確地得到.已有許多數(shù)值方法被提出用于對高維NLS/GP方程進行模擬,比如分裂格式有限差分法、有限元法、時間偽譜法、修正歐拉算法、蒙特卡羅法和無網格方法等[3-14],然而上述方法屬于網格類或是基于背景網格,對非均勻節(jié)點布置下或高維復雜區(qū)域問題的模擬實現(xiàn)都較復雜[15-18].為此已有許多學者對完全不依賴于網格的粒子方法(或純無網格方法)進行了關注,并成功將粒子方法推廣應用到實數(shù)域上偏微分方程的求解,如光滑粒子動力學(smoothed particle hydrodynamics,SPH)方法[15-25]等.然而上述純無網格方法對非線性薛定諤方程的研究還處于起步階段,特別對高維復雜GP方程的研究在國際上尚不多見.

作為完全不依賴于網格的SPH粒子方法,在高維非線性薛定諤方程上的模擬應用還處在初步研究階段[15],其原因在于傳統(tǒng)SPH方法的數(shù)值精度低、穩(wěn)定性差[16-19,20-22].SPH方法精度和穩(wěn)定性的改進與高維非線性薛定諤方程的數(shù)值研究是國際上的兩個熱點問題.針對SPH方法數(shù)值精度的提高,已有一些改進SPH方法[16-19,21-25]被提出,并被應用于許多領域[16-25],但它們仍存在一些自身缺陷,因此當推廣應用到非線性薛定諤方程的模擬時需要進一步改進.粒子方法模擬高維NLS/GP方程較網格類方法的主要優(yōu)勢在于：不依賴于網格可以任意布點;便于處理復雜區(qū)域和計算空間導數(shù);模擬程序易于實現(xiàn),并且便于實施并行計算以提高計算效率.

基于上述分析,本文針對含有非線性項和旋轉導數(shù)項高維NLS/GP方程的準確、高效性粒子法模擬給出一種四階分裂修正并行SPH (HSSCPSPH)方法.提出的HSS-CPSPH方法思想主要來源于：首先,采用分裂思想將NLS方程分解為線性導數(shù)項和非線性項兩個微分方程;其次,基于Taylor展開拓展應用文獻[17,21]中修正SPH方法對含線性導數(shù)方程進行二階精度顯式求解;再次,運用四階對稱分裂格式對上述兩個方程進行積分離散;最后引入基于粒子搜索的MPI并行技術提高計算效率.給出的HSS-CPSPH方法將兼具傳統(tǒng)SPH方法和已有改進SPH以及分裂格式的優(yōu)點,也較已有分裂有限差分法[6-8]具有更好、更靈活的模擬應用性.在數(shù)值算例模擬中,將HSSCPSPH方法的數(shù)值結果與解析解或其他數(shù)值結果進行了比較,分析了數(shù)值精度、收斂階和計算效率.比較分析結果表明,給出的HSS-CPSPH方法可以準確、高效地求解周期和Dirichlet邊界條件下二維或三維NLS/GP方程,并能成功地預測二分量孤立波傳播中的非線性奇異特性和帶外旋轉算子下量子渦旋的復雜變化過程.

2 非線性薛定諤/Gross-Pitaevskii方程

量子力學中許多非線性物理現(xiàn)象常用Gross-Pitaevskii (GP)方程[3,4,6,7]來描述,其高維直角坐標系下無量綱化的控制方程為

γy,γz是實常數(shù);Lz=-i(x?y-y?x) 是玻色-愛因斯坦凝聚態(tài)(Bose-Einstein condensates,BEC)下無量綱速度旋轉項Ω的角動量算符的z分量.

為封閉上述方程,考慮如下初始條件：

邊界條件

或周期邊界條件φ(x,t)=φ(x+l,t),其中l(wèi)表示方程周期,或φ(x,y,t)=φ(x+l1,y,t),φ(x,y,t)=φ(x,y+l2,t),其中 (l1,l2)為 (x,y) 方向的周期.

根據(jù)文獻[6—8]可知,上述無窮區(qū)域邊界條件(4)式可以近似處理為齊次Dirichlet邊界,周期邊界條件采用兩邊施加幾層虛粒子方式進行處理(詳見3.3節(jié)).

3 高階分裂修正并行SPH方法

針對高維非線性GP方程的粒子方法求解,拓展應用文獻[17,21]中給出的一階修正SPH方法時,隨著模擬時間的延長會出現(xiàn)精度低和穩(wěn)定性差的現(xiàn)象,這是因為本文考慮的GP方程中存在非線性項和一階導數(shù)項,對模擬算法的精度和穩(wěn)定性要求較高.為此,在時間高階分裂法[1,2]和已有改進SPH方法[15-18]的基礎上,引入MPI并行算法,提高高維問題模擬計算效率.本文對二維/三維GP方程模擬給出一種HSS-CPSPH方法,基本思想是：首先,采用時間分裂思想將GP方程分成線性導數(shù)和非線性項兩個方程;其次,對線性導數(shù)方程拓展應用文獻[15,17]中的改進SPH方法進行離散;再次,引入文獻[1,2]中四階分裂格式對分裂的兩個方程進行交替求解;最后,采用基于MPI的并行算法來提高計算機模擬效率.

3.1 時間分裂法[1-3,6]

作為一種精確高效的算法,時間分裂法在復數(shù)域非線性薛定諤方程求解中已被廣泛應用,根據(jù)文獻[6—8],GP方程(1)可以重寫為

和非線性問題

上述過程的成立條件和精度可見文獻[26].

本文對方程(6)和(7)的求解采用如下四階分裂法[1,2]：第一步求解方程(7),第二步用第一步得到的解作為初始條件來求解方程(6),第三步用第二步得到的解作為初始條件求解方程(7),第四步使用第三步得到的解作為初始條件求解方程(6),第五步用第四步得到的解作為初始條件求解方程(7),此時得到的解作為下一時間層的數(shù)值解.

3.2 高階分裂修正SPH離散格式

3.2.1 修正SPH離散格式

在傳統(tǒng)SPH方法模擬中,首先使用核函數(shù)進行積分插值,其次在模擬區(qū)域Ω內,對有限個粒子賦予密度、質量、溫度等物理含義,最后用粒子離散化方式近似得到核函數(shù)的積分形式.

在任意點x= (x,y,z)處,函數(shù)φ(x) 及其微分形式可以用核函數(shù)W和核函數(shù)的導數(shù)?W表示[19]：

其中h為光滑長度,決定了核函數(shù)W的支持域范圍,W一般要求滿足正則化、對稱性、緊致性和Dirac函數(shù)性質[19].本文模擬中取h=0.9d0-1.1d0(d0為初始粒子間距),取分段五次樣條核函數(shù)[19],對應支持域尺寸為 3h.

引入體積v=m/ρ(其中m,ρ都是實常數(shù)),對于任意位置xi=(xi,yi,zi)處粒子i,其支持域內的相鄰粒子為j,則有粒子近似公式：

Wij=W(|xi-xj|,h) ,?iWij=?W(|xi-xj|,h)/?xi,且?iWij滿足關系

為了提高傳統(tǒng)SPH方法的精度和穩(wěn)定性,基于Taylor展開,人們已經提出了一些改進SPH方法[15-22].本文拓展應用文獻[17,21]中具有二階精度的一階修正對稱SPH方法,該方法首先將二階導數(shù)項分為兩個一階導數(shù),然后對一階導數(shù)基于Taylor展開進行離散,則(11)式變?yōu)?/p>

3.2.2 高階分裂修正SPH格式

將3.1節(jié)四階分裂格式與3.2.1節(jié)修正SPH格式進行結合,對方程(1)進行離散,可得如下高階分裂修正SPH離散格式

其中i= 0,1,2,···,M;n=0,1,···;M是區(qū)間內所有粒子總數(shù),n是時間層,dt是時間步長,

四階分裂格式參數(shù)

上述離散格式被稱為高階分裂修正SPH方法,結合3.4節(jié)MPI并行計算技術將本文提出的方法稱為“HSS-CPSPH”,其中對線性導數(shù)方程(6)的時間離散采用二階龍格-庫塔顯格式,根據(jù)文獻[16—19]中的穩(wěn)定性限制條件,時間步長取 dt≤.本文方法與文獻[15]中給出的SS-ICPSPH方法主要的不同之處在于：對線性導數(shù)方程的離散本文采用二階龍格-庫塔顯格式,后者采用隱格式;本文采用四階分裂格式,后者采用二階分裂格式,使得本文方法精度較文獻[15]中的方法精度高(見第4節(jié));文獻[15]運用提出的方法僅對帶Dirichlet邊界的GP方程進行模擬研究,本文運用給出的方法對帶周期邊界和Dirichlet邊界的非線性薛定諤問題都進行了模擬分析,并且準確預測了周期邊界下孤立波傳播中出現(xiàn)的奇異現(xiàn)象(見5.1節(jié)).

3.3 邊界條件處理

運用上述提出的方法對高維GP方程進行模擬時,初邊值條件的準確處理對數(shù)值模擬的精度和穩(wěn)定性至關重要.對初始條件可以準確離散施加,邊界條件(4)式可以近似處理為Dirichlet邊界.以二維區(qū)域情況為例,周期性邊界條件處理如下：

其中Xa,Xb,Ya,Yb分別為x,y方向區(qū)域的最小和最大值.將區(qū)域沿 (x,y)方向分別均分成J和K整數(shù)份,令空間步長h1=(Xb-Xa)/J,h2=(Yb-Ya)/K,xj=Xa+jh1,0≤j≤J-1 ,yk=Ya+kh2,0≤k≤K-1 .為防止粒子方法處理周期邊界的粒子缺失問題,在區(qū)間外側各擴充4個點(大于等于 3h),可以得到

則具有 (l1,l2) 周期邊界條件的方程滿足如下條件：

本文在第4和第5節(jié)通過算例模擬驗證了上述邊界處理的準確性和有效性.

3.4 粒子搜索并行

SPH粒子方法模擬中,相鄰粒子搜索和高維問題下需要幾十萬至上百萬粒子物理量更新計算,會使得計算機模擬中占較大計算內存和較長CPU計算時間,因此運用提出方法對高維GP進行模擬時提高計算效率是必要的.根據(jù)模擬中粒子位置不變的特點,引入文獻[17,21]中基于MPI粒子搜索并行技術.即基于MPI并行算法,首先對每個粒子支持域內的相鄰粒子進行標記,標記后對所有粒子進行物理量賦值,主要是質量和密度,為了提高計算效率,在并行計算中,不同節(jié)點間會進行密度通信.在計算過程中節(jié)點間進行空間上的一階導數(shù)項通信以及最后計算函數(shù)值的更新.本文的并行計算是將所有粒子進行編號,再根據(jù)CPU數(shù)量將粒子平均分配,進行同時計算.每個粒子與之相互作用的相鄰粒子支持域范圍通過下列核函數(shù)確定,

其中q=rij/h,rij=|xi-xj|,系數(shù)w0為 120/h(一維),7/(478πh2)(二維),3/(359πh3)(三維).本文中的支持域范圍是以 3h為半徑的圓或球.

4 HSS-CPSPH方法的精度和收斂性分析

本節(jié)通過對具有解析解帶不同邊界條件的非線性薛定諤方程進行模擬,對提出的HSSCPSPH方法的精度、收斂階以及計算效率進行分析.分別對二維帶周期邊值、一維二分量帶周期邊值和二維帶第一類邊值的問題進行數(shù)值模擬,并與解析解做比較,分析所提方法的準確性和高效性,與已有SS-ICPSPH方法[10]做對比證明本文所提方法具有較高精度和較快的收斂速度.

為分析方法的數(shù)值精度和收斂速度,定義如下最大誤差范數(shù)和收斂階：

其中d01和d02代表了兩種不同的粒子初始間距.

4.1 二維帶周期邊界下孤立波傳播

考慮方形區(qū)域 [0,2π]×[0,2π] 上帶周期邊界的非線性薛定諤方程[3]

對應解析解為

u(x,y,t)=u(x+2π,y,t),

u(x,y,t)=u(x,y+2π,t),

(x,y)∈R×R,0＜t≤T,T=1.

該算例模擬中,k1=k2=1時采用 129×129個均勻分布粒子,時間步長為 dt=10-4(見圖1);k1=k2=4時采用h= π/32 個均勻分布粒子,時間步長為 dt=10-4(見圖2).圖1和圖2分別展示了k1和k2取不同系數(shù)時由HSS-CPSPH方法獲得的u(x,π) 的實部曲線圖,并將其與解析解及SSICPSPH結果進行對比,可以看出,HSS-CPSPH結果與解析解相符合,但隨著時間的延長,兩種數(shù)值結果產生較小的偏差,這與文獻[3]中分裂有限差分法得到數(shù)值模擬結論類似.

為體現(xiàn)提出的HSS-CPSPH方法求解周期性問題的精度和收斂速度,表1和表2分別列出了模擬較短時間內數(shù)值結果的誤差和收斂階,由表1和表2可以看出：1)em隨著粒子數(shù)的增加而減小,HSS-CPSPH方法得到的誤差較SS-ICPSPH方法的小;2) HSS-CPSPH方法較SS-ICPSPH方法具有更快的收斂速度.為進一步體現(xiàn)粒子方法在粒子分布非均勻情況下數(shù)值模擬的精度,表3列出了粒子分布均勻和兩種非均勻情況下(見文獻[15])兩個不同時刻的最大誤差em.由表3可知,粒子方法在粒子分布均勻和分布非均勻情況下得到的數(shù)值結果誤差都比較接近,表明了該方法易推廣應用到非規(guī)則區(qū)域問題的模擬,且保持較好的精度.

表1 k1=k2=1,h=π/64 時幾個不同時刻里兩種方法的誤差emTable 1.Erroremobtained using two different methods at different time (k1=k2=1,h=π/64 ).

圖1 k1=k2=1,h=π/64 時不同時刻u(x,π) 的實部沿x軸的變化 (a)t = 1;(b)t=3Fig.1.Curve of the R e(u(x,π)) alongx-axis at different time withk1=k2=1,h=π/64 ：(a)t = 1;(b)t = 3 .

圖2 k1=k2=4,h=π/128 時不同時刻u(x,π) 的實部沿x軸的變化 (a)t = 0.1;(b)t = 1Fig.2.Curve of the R e(u(x,π)) alongx-axis at different time withk1=k2=4,h=π/128 ：(a)t = 0.1;(b)t = 1.

表2 k1=k2=1 ,時間t = 1時,兩種方法在不同粒子間距下的誤差和收斂階Table 2.Erroremand convergent orderorαobtained using two different methods att = 1 and different particle distance (k1=k2=1 ).

表3 k1=k2=1,h=π/64 時,粒子分布均勻或不均勻方式下,兩種方法的誤差emTable 3.Erroremobtained using different methods at different distribution (k1=k2=1 ,h = π/64 ).

4.2 一維二分量帶周期邊界下孤立波傳播

為進一步展示本文提出的HSS-CPSPH方法能夠準確捕捉多分量多孤立波傳播中的尖角現(xiàn)象,本小節(jié)考慮區(qū)間 [-20,80] 上一維二分量帶周期邊界兩種初始條件下的非線性薛定諤方程[14]

初始條件1為

其中參數(shù)α= 0.5,β=2/3,c=1,a=1,x0=0 .

具有三個孤子波傳播的初始條件2為

其中α= 0.5,β= 2/3,c1=1,c2=0.1,c3=-1 ,a1=1,a2=0.72,a3=0.36,x1=0,x2=25 ,x3=50 .

周期邊界條件u(x,t)=u(x+100,t) ,上述兩種初始條件下對應的解析解可參見文獻[14].

圖3和圖4給出了兩種不同初始條件下幾個時刻里分量u的不同數(shù)值結果.觀察圖3和圖4可知,HSS-CPSPH結果與解析解相符合,即使在較長模擬時間里;提出的HSS-CPSPH方法能準確捕捉到孤立子傳播中奇異現(xiàn)象;HSS-CPSPH方法可以準確得到一定時間里孤立子波的傳播過程.

4.3 二維第一類邊界下孤立波傳播

為驗證提出的HSS-CPSPH方法求解帶第一類邊界非線性薛定諤問題的準確性,以及進一步體現(xiàn)它與文獻[15]給出的粒子方法相比具有的優(yōu)點.本小節(jié)考慮區(qū)域 [0,2π]2內的二維非線性薛定諤方程[8]

初值條件為u0(x,y)=sinxsiny,第一類邊界條件可以通過方程解析解得到,其中V(x,y)=1-sin2xsin2y.對應方程解析解為uexact(x,y,t)=sinxsinyexp(-i2t) .

圖3 初始條件1下,在4個不同時刻孤立波函數(shù)|u|的傳播過程 (a)t = 0;(b)t = 20;(c)t = 30;(d)t=50Fig.3.Solitary wave propagation process of|u|at different time with the initial condition 1：(a)t = 0;(b)t = 20;(c)t = 30 ;(d)t = 50 .

圖4 初始條件2下,在4個不同時刻三孤立子波函數(shù)|u|的傳播過程 (a)t = 0;(b)t = 20;(c)t = 30;(d)t=50Fig.4.Solitary wave propagation process of|u|at different time with initial condition 2：(a)t = 0;(b)t = 20;(c)t = 30 ;(d)t = 50 .

表4 h = π/64 時,三個不同時刻兩種方法的最大誤差emTable 4.Erroremobtained using two different methods at three times (h = π/64 ).

表4和表5列出了兩種粒子方法模擬上述方程的誤差和收斂速度.通過表4和表5可知,在模擬第一類邊值NLS方程時,本文提出的HSSCPSPH方法較文獻[15]中的粒子方法具有較小誤差和較快收斂速度.此外,為了體現(xiàn)本文算法通過并行提高計算效率的必要性,模擬過程中基于MPI并行算法,采用了多個CPU進行計算機模擬,取粒子數(shù)為66049( 2 572),時間步長為 dt=5×10-5.通過實際計算模擬知,用一個CPU串行計算時,運行第1步的時間約為49.89 s,之后100步平均每步約1.5 s;當采用4個CPU時,運行第1步的時間約為12.0586 s,之后100步平均每步約0.25 s.可以看出,運行第1步時4個CPU的計算時間基本上是1個CPU的計算時間的1/4,計算步數(shù)增加后,4個CPU平均每步的計算時間略小于1個CPU平均每步計算時間的1/4.這與理想的1/4有一定偏差,可能原因有下面幾點：1) 同一節(jié)點上不同CPU之間通訊消耗時間非常小;2) 此處為二維模擬區(qū)域粒子數(shù)不是很多的情況;3) 得到的計算時間有舍入誤差.當模擬粒子數(shù)增大(如5.2節(jié)三維區(qū)域模擬中涉及幾百萬及以上粒子),CPU個數(shù)增加到在不同節(jié)點上時,會使得平均每步的計算時間增加,不同節(jié)點上通訊消耗時間延長,此時經實際模擬得到的結論會與實際計算理論相符合,即隨著CPU數(shù)量增加,計算效率的提高倍數(shù)要明顯小于CPU增加的倍數(shù)(詳見5.2節(jié)并行計算效率分析).

表5 t = 1時不同空間步長情況下兩種粒子方法的誤差和收斂階Table 5.Error and order of convergence by different methods att = 1 and differenth.

5 GP方程數(shù)值模擬

為了更好地體現(xiàn)提出的基于分裂格式的粒子方法模擬NLS/GP問題的高效、準確性,本節(jié)選取了一個二維周期性無解析解NLS問題和兩個BEC下帶外旋轉項的GP問題,分別預測了周期邊界下孤立子傳播的奇異特性和BEC下量子化渦旋變化過程,并將模擬結果與其他數(shù)值方法[8,15]計算結果做對比.模擬中,也通過三維算例的計算機模擬展示了采用基于MPI并行算法提高計算效率的必要性.

5.1 二維周期邊界無解析解NLS方程

考慮區(qū)間 [0,2π]×[0,2π] 上的非線性薛定諤方程[3]

初值條件為u(x,y,0)=(1+sinx)(2+siny) ,參數(shù)β=1,周期邊界長度為 (2π,2π) .

上述NLS方程是帶初始和周期邊界下孤立子傳播中出現(xiàn)奇異特性的典型問題,目前該方程無解析解,常被用來驗證提出算法模擬帶周期邊界下孤立子傳播過程的可靠性(見文獻[3]).圖5給出了初始時刻和t= 0.108 時刻三種數(shù)值方法的模擬結果.從圖5可知,在t= 0.108 時出現(xiàn)明顯的孤立子波奇異特性,三種數(shù)值結果相符合,進一步表明本文提出的粒子方法捕捉孤立子波傳播中奇異現(xiàn)象是可靠的.

5.2 三維Gross-Pitaevskii方程

考慮三維BEC下帶旋轉項的Gross-Pitaevskii方程[15]

初始條件為

其中,V(x,y,z,t)=(γxx2+γyy2+γzz2)/2 ,參數(shù)Ω=0.75,(γx,γy,γz,β) = (1.2,0.8,2,50) .

為表明HSS-CPSPH方法模擬三維GP問題的可靠性和展示BEC下不同截面上波函數(shù)值隨時間的變化過程,圖6給出了沿y軸方向上兩種不同數(shù)值曲線,圖7給出了由HSS-CPSPH方法得到的3個不同時刻沿u(0,y,z) 和u(x,y,0) 兩個截面上波函數(shù)變化等值線.圖6中本文方法得到的變化曲線與HO-SSFDM方法[15]結果一致,體現(xiàn)了HSS-CPSPH方法模擬預測三維GP問題是有效的;波函數(shù)的峰值隨時間演化出現(xiàn)先變小后變大的趨勢,但波峰值的總體變化趨勢是慢慢變小.通過觀察圖7可知,不同截面上隨時間演化波函數(shù)值的分布是不同的、復雜的,且等值線峰值隨時間延長逐漸變小.

圖5 在兩個不同時刻不同數(shù)值方法得到的|u|等值線圖(a)t = 0;(b)t=0.108Fig.5.Contours of |u| obtained using different methods at two different times：(a)t = 0;(b)t = 0.108 .

圖6 不同時刻|u|2沿y軸 (x=0,z=0) 變化曲線Fig.6.Curve of |u|2 alongy-axis (x=0,z=0) at different time.

圖7 在3個不同時刻|u|2在不同截面上的等值線 (a)(0,y,z) 截面;(b) (x,y,0) 截面Fig.7.Contour of |u|2 along different profile at different time：(a)(0,y,z);(b)(x,y,0) .

為進一步體現(xiàn)運用本文粒子方法模擬高維GP問題時引入MPI并行算法以減少CPU消耗時間的必要性,對固定粒子數(shù)和增加粒子數(shù)兩種情況采用不同CPU個數(shù)的計算效率進行了分析(見表6和表7).本文并行計算機采用Red Hat Enterprise Linux 5.8 × 86_64操作系統(tǒng),MPI類型為Intel MPI Toolkits 4.0.3,總共有17個IBM BladeCenter HS22計算節(jié)點,每個節(jié)點包括12個兩路六核CPU,CPU型號為Intel Xeon 6C X5650,主頻為2.67 GHz.表6列出了固定粒子數(shù) 1613(約四百多萬)增加CPU個數(shù)情況下運行不同步數(shù)時所消耗的CPU時間,表7列出了不同粒子數(shù)下不同CPU情況下,運行到1000步時平均每步所消耗的CPU時間.為了更加直觀地看出并行計算對計算效率的提高,此處定義相對加速比S為單節(jié)點上并行算法的運行時間/N個節(jié)點上并行算法的運行時間(每個節(jié)點上12個CPU).選用表6中12,24,36,72個CPU在計算到1000步時的時間作為研究對象,分別對應著1,2,3,6個節(jié)點.

由表6和表7可以看到,隨著CPU個數(shù)的增加,計算效率提高很快,但計算效率提高的比值是低于CPU個數(shù)增加的比值,這是因為CPU個數(shù)增加后不同CPU上得到結果相互通訊的時間也將增加.從加速比值上可以明顯看出,在開始使用1,2,3個節(jié)點時,計算效率的提升是十分明顯的,但略低于CPU增加的比值,這是由于不同節(jié)點之間通訊消耗一定的時間導致的,但是當節(jié)點數(shù)繼續(xù)增加到6個節(jié)點時,通訊消耗時間也隨之增加,使得計算效率增加減緩,明顯低于CPU增加的比值,這與計算機相關計算理論也較為符合;計算機模擬中運行第一步消耗CPU時間比后續(xù)運行的平均步消耗CPU時間要長,且隨粒子數(shù)增加這個消耗CPU時間的比值要變大,這是因在計算物理量循環(huán)更新前需要標定每個粒子的相鄰粒子,此標定時間隨粒子數(shù)增加而增加;在CPU個數(shù)不變的情況下,隨著粒子數(shù)增加平均每個計算步消耗CPU時間也將變長,且該變長時間比值與粒子數(shù)增加比值是呈非線性的.由于本文采用的并行計算主要是通過增加CPU個數(shù)來提高計算效率,不同CPU數(shù)下的數(shù)值近似算法不變,使得多個CPU下的計算結果與單個CPU的計算結果是一致的,通過不同CPU下的實際模擬結果對比發(fā)現(xiàn)亦是如此,從而表明多CPU下的求解質量是可靠的.

表6 粒子數(shù)為 1613 時,不同CPU個數(shù)下運行到不同步數(shù)所需時間(單位：s)Table 6.Consumed CPU time (unit：s) of different calculated time step with particle number 1613 at different CPUs.

表7 在不同粒子數(shù)下不同CPU個數(shù)下,運行到1000步時平均每步所消耗時間(單位：s)Table 7.The average consumed CPU time (unit：s)of calculated time step 1000 with different particle number and different CPUs.

通過上述模擬分析得知,本文基于MPI并行計算給出的HSS-CPSPH方法模擬三維GP問題是高效、可靠的.

5.3 二維二分量Gross-Pitaevskii方程

為進一步展示提出的HSS-CPSPH方法預測BEC下量子渦旋隨時間演化過程的可靠性,本小節(jié)考慮區(qū)間 [-8,8]2上二維二分量BEC下帶旋轉項的Gross-Pitaevskii方程[8]

初邊值條件為

圖8 兩個不同時刻下|u1|沿x軸(y = 0.5)的變化 (a)t = 0.05;(b)t = 0.25Fig.8.Curve of |u1| alongx-axis (y = 0.5) at two different time：(a)t = 0.05;(b)t = 0.25.

圖9 兩 個不同時刻 下t = 0 (第一列)和t = 0.25 (第二列)三個物理量等值線變化 (a1),(a2) Re(u1);(b1),(b2) I m(u1) ;(c1),(c2)|u1|Fig.9.Contours of three physical quantities at two different timest = 0 (the first row) andt = 0.25 (the second row)：(a1),(a2) Re(u1) ;(b1),(b2) I m(u);(c1),(c2)|u|.

其中,wi(x,y)=1.5x2+0.5y2,i=1,2,Θ=0.7,σ=-100,?=0.8 .

圖8給出了兩個時刻兩種不同數(shù)值方法得到的波函數(shù)沿x軸(y= 0.5)的變化,圖9給出了由HSS-CPSPH方法得到的兩個時刻3個物理量(第一分量u1的實部、虛部和模)等值線.由圖8和圖9知,HSS-CPSPH結果與SS-FDM結果[7]相符,表明本文粒子法模擬二分量GP問題是可靠的;HSS-CPSPH方法能夠準確預測出具有二分量BEC下帶旋轉項量子渦旋隨時間的變化過程,進一步展示了所提方法高效、準確預測GP問題的能力.

6 結 論

本文為提高SPH方法模擬含孤立波非線性問題的數(shù)值精度和穩(wěn)定性,將四階精度時間分裂格式與修正SPH方法相結合,并引入基于粒子搜索的MPI并行計算技術,給出一種能夠高效準確地模擬高維NLS/GP方程的HSS-CPSPH方法.數(shù)值模擬中,考慮了粒子分布均勻和非均勻情況下帶有兩種不同邊界條件的二維和三維NLS/GP方程,并與解析解或其他數(shù)值結果做對比,分析了所提方法的數(shù)值精度、收斂階及數(shù)值預測的可靠性.所有數(shù)值算例表明：

1)提出的HSS-CPSPH方法求解高維非線性薛定諤方程較已有粒子方法具有較高精度和較快收斂速度,且較網格類方法具有更好、更靈活的應用性;

2)無論在粒子分布均勻還是非均勻情況下,HSS-CPSPH方法模擬NLS方程都具有較高的數(shù)值精度;

3)給出的HSS-CPSPH方法成功地預測了周期邊界下二維NLS方程描述的孤立波傳播中的奇異特性,準確地展示了BEC下帶外旋轉項三維單分量和二維二分量GP方程中量子渦旋隨時間演化過程.

值得注意的是,目前未見文獻將四階分裂格式與修正SPH方法耦合,并引入MPI并行技術對NLS/GP方程進行模擬研究.本文對不同邊界下高維NLS方程的模擬給出的HSS-CPSPH法較已有粒子方法具有較高精度和較快收斂速度,較網格類方法具有更好、更靈活的拓展應用性,為高維NLS方程的數(shù)值模擬提供了一種高效準確的粒子方法.