岳孝強,　周志陽,　徐小文,　舒　適

(1.湘潭大學(xué) 數(shù)學(xué)與計算科學(xué)學(xué)院,湖南 湘潭 411105;2.北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100094)

在慣性約束聚變(ICF)等高能量密度物理領(lǐng)域的數(shù)值模擬中，輻射流體力學(xué)方程組的求解是非常重要的環(huán)節(jié).該方程組由描述流體運動與輻射傳輸?shù)膬深惙匠恬詈隙?，其中后者在擴散近似下由三溫輻射擴散方程組(RDEs)刻畫[1].由于物理量間的復(fù)雜耦合等特點，三溫RDEs離散系統(tǒng)的條件數(shù)極差，致使求解時間通常占整體模擬時間的70%以上[2-3].因此，研究三溫RDEs的高效并行解法器是一項具有重要理論意義與實用價值且難度很大的課題.預(yù)條件Krylov子空間迭代法是求解三溫RDEs離散系統(tǒng)的主要方法，其中影響算法效率的關(guān)鍵是預(yù)條件子的構(gòu)造，目前已有很多研究工作[4-10].值得一提的是，[9][10]針對Euler網(wǎng)格構(gòu)造了基于物理量粗化及耦合強度的PCTL法及其自適應(yīng)方法(但限于串行實現(xiàn)).這些預(yù)條件算法(主要是ILU類、多重網(wǎng)格類以及兩者的有機組合)和理論主要是針對二維結(jié)構(gòu)和協(xié)調(diào)非結(jié)構(gòu)網(wǎng)格，仍需不斷完善與發(fā)展，特別是關(guān)于算法的并行設(shè)計以及可擴展性仍需進(jìn)一步提升.

自適應(yīng)結(jié)構(gòu)網(wǎng)格應(yīng)用支撐框架(JASMIN)以其數(shù)據(jù)結(jié)構(gòu)、高效算法以及并行通信等高度封裝的特性而成為大規(guī)模實際應(yīng)用需求的一個重要研發(fā)平臺[11].基于細(xì)化模式的隱式時間積分算法是求解ICF內(nèi)爆過程中多物理耦合問題的常用模式[12].在該模式下，求解網(wǎng)格片層次結(jié)構(gòu)上離散系統(tǒng)的主要工作量體現(xiàn)在每個網(wǎng)格層上離散系統(tǒng)的求解.因此，如何在不含懸點的分片結(jié)構(gòu)網(wǎng)格上，為三溫RDEs的大規(guī)模離散系統(tǒng)設(shè)計出基于自適應(yīng)PCTL預(yù)條件子的高效并行PGMRES解法器，將為解決基于細(xì)化模式的ICF瓶頸問題起到積極作用.

本文針對三溫輻射擴散問題的全隱有限體積格式，將自適應(yīng)PCTL法推廣到三維情形，給出已有的二維自適應(yīng)PCTL法在JASMIN平臺下的并行實現(xiàn)算法(涉及進(jìn)程分組及數(shù)據(jù)轉(zhuǎn)換)，并通過二維和三維典型算例驗證了新并行解法器的高效性與良好的算法可擴展性.

1　模型問題及有限體積格式

考慮三溫RDEs

(1)

(1)是非線性偏微分方程組，利用向后Euler方法進(jìn)行時間離散、“凝固系數(shù)法”線性化，再采用全隱有限體積格式，可得大規(guī)模離散系統(tǒng)Au=f，其中A∈R3n×3n，n為網(wǎng)格規(guī)模.假定自由度按照光子、電子、離子進(jìn)行排列，則A可寫為塊三對角結(jié)構(gòu)

(2)

2　一種基于JASMIN平臺的并行自適應(yīng)PCTL法

2.1　一種自適應(yīng)PCTL法

(3)

需要注意的是，自適應(yīng)PCTL法的粗化過程是基于物理量進(jìn)行的，算法描述.

算法1 自適應(yīng)PCTL法:w=Bαg對給定的兩組參數(shù)θwc、σwc、θsd、σsd;ε、μc、να、μα(α=R,E,I),若滿足|{i:-dERii≤θwc×aEii,i=1,2,…,n}|/n≥σwc,(4)|{i:-dEIii≤θwc×aEii,i=1,2,…,n}|/n≥σwc,(5)則w=A～-1R000A～-1E000A～-1Ié?êêêêêù?úúúúúg其中A～-1α(α=R,E,I)為A-1α的近似,具體做法:若|{i:jaαij≥θsd×aαii,i=1,2,…,n}|/n≥σsd,(6)則調(diào)用να次Jacobi迭代求解;否則采用AMG法求解.若不滿足式(4)或式(5),則作以下三步:前磨光過程:分別求解電子、輻射和離子方程子系統(tǒng)AEwE=gE,ARwR=gR-DREwE,AIwI=gI-DIEwE.具體做法:若Aα滿足式(6),則采用Jacobi迭代求解;否則采用AMG法求解,取最大迭代次數(shù)、迭代控制精度分別為μα、ε.求解粗水平子系統(tǒng):(PTAP)wc=PT(g-Aw),其中P由式(3)計算.具體做法:采用AMG法求解,取最大迭代次數(shù)、迭代控制精度分別為μc、ε.粗水平校正:w=w+Pwc.

2.2　自適應(yīng)PCTL法的并行實現(xiàn)

由算法1可知：自適應(yīng)PCTL預(yù)條件子涉及若干子系統(tǒng)的求解，且其中有些子系統(tǒng)可以獨立求解.基于這個發(fā)現(xiàn)，本節(jié)將利用進(jìn)程分組策略給出算法1的并行實現(xiàn).為此，首先給出基于進(jìn)程分組的分劃數(shù)組的概念.

(7)

為下標(biāo)相匹配，以下用記號GR、GE、GI來分別表示G1、G2、G3.

在JASMIN中，式(2)的并行數(shù)據(jù)為各子塊的并行矩陣和交叉塊的并行向量：

AR(pu),AE(pu),AI(pu)；VRE(pu),VER(pu),VEI(pu),VIE(pu)，

(8)

為便于算法1基于進(jìn)程分組策略的實現(xiàn)，需要對數(shù)據(jù)式(8)進(jìn)行轉(zhuǎn)換，此時需指定分劃數(shù)組p，這里我們將其取為p(j)=pu(3j),j=0,1,…,m.數(shù)據(jù)轉(zhuǎn)換的目標(biāo)是生成

(1) 并行數(shù)據(jù)I：各子進(jìn)程組上的并行數(shù)據(jù)

GR:AR(p),VRE(p)；GE:AE(p),VER(p),VEI(p)；GI:AI(p),VIE(p)，

式中：Aα(p)(α=R,E,I)表示子矩陣Aα在子進(jìn)程組Gα上基于分劃數(shù)組p的并行矩陣；Vαβ(p)(α,β=R,E,I)表示Vαβ在子進(jìn)程組Gα上基于分劃數(shù)組p的并行向量.

(2) 并行數(shù)據(jù)II：進(jìn)程組G上基于分劃數(shù)組pS的并行矩陣A(pS).

上述并行數(shù)據(jù)轉(zhuǎn)換的本質(zhì)是將矩陣和向量數(shù)據(jù)在各進(jìn)程上進(jìn)行分配，其關(guān)鍵是如何建立接口并行數(shù)據(jù)和目標(biāo)并行數(shù)據(jù)之間的映射關(guān)系.以下分兩步來完成上述并行數(shù)據(jù)轉(zhuǎn)換.

① 將并行數(shù)據(jù)式(8)轉(zhuǎn)換為并行數(shù)據(jù)I.(a) 將并行矩陣Aα(Pu)轉(zhuǎn)為Aα(p),α=R,E,I.具體的數(shù)據(jù)分配方式為：將Aα(pu,3i)，Aα(pu,3i+1)，Aα(pu,3i+2)分配給進(jìn)程Gα(i)，這里Aα(pu,k)表示并行矩陣Aα(pu)屬于G(k)的長方陣.(b) 將并行向量Vαβ(pu)轉(zhuǎn)為Vαβ(p)，α,β=R,E,I.具體的數(shù)據(jù)分配方式為：將Vαβ(pu,3i)，Vαβ(pu,3i+1)，Vαβ(pu,3i+2)分配給進(jìn)程Gα(i)，這里Vαβ(pu,k)表示并行向量Vαβ(pu)屬于G(k)的子向量.

算法2 并行自適應(yīng)PCTL法的Setup過程Step1 將進(jìn)程組G分成子進(jìn)程組GR、GE、GI.Step2 生成弱耦合條件和強對角占優(yōu)條件是否成立的標(biāo)志變量.具體為子進(jìn)程組Gα(α=R,I):利用Aα(p),按式(6),生成強對角占優(yōu)標(biāo)志標(biāo)量ddα;子進(jìn)程組GE:利用AE(p)和VER(p),VEI(p),按式(4)至式(6),生成強對角占優(yōu)標(biāo)志變量ddE和弱耦合條件標(biāo)志變量wcE,并將wcE廣播給其他兩個子進(jìn)程組.Step3 若wcE=1,則子進(jìn)程組Gα(α=R,E,I)分別作以下步驟:若ddα=0,則對Aα(p)調(diào)用BoomerAMG法[13]的Setup算法.Step4 若wcE=0,則執(zhí)行以下步驟4.1 三個子進(jìn)程組分別作以下步驟(1)子進(jìn)程組Gα(α=R,I):①對Aα(p)調(diào)用BoomerAMG法的Setup算法;②調(diào)用BoomerAMG法的Solve算法求解:Aα(p)pα(p)=-VαE(p),得pα(p).(2)子進(jìn)程組GE:若ddE=0,則對AE(p)調(diào)用BoomerAMG法的Setup算法.4.2 利用pR(p),pI(p)和分劃pu生成并行插值矩陣P(pS).4.3 利用Aα(pu),VRE(pu),VRE(pu),VER(pu),VEI(pu)以及pR(p),pI(p)生成并行粗化矩陣Ac(pu).4.4 對Ac(pu)調(diào)用BoomerAMG法的Setup算法.

算法3 并行自適應(yīng)PCTL法的Solution過程:w=Bαg若wcE=1,則子進(jìn)程組Gα(α=R,I)分別作以下步驟,生成w(pS):以0為初值,對Aα(p)wα(p)=gα(p)作να次迭代,其中:當(dāng)ddα=1時采用并行Jacobi法作迭代,否則采用BoomerAMG法的Solve算法作迭代.若wcE=0,則執(zhí)行以下步驟:調(diào)用ParaRelaxBlockGS算法,作一次塊Gauss-Seidel磨光,生成w(pS);計算殘量:r(pS)=g(pS)-A(pS)w(pS);殘量限制:rc(pu)=P(pS)Tr(pS);調(diào)用BoomerAMG法求解Ac(pu)wc(pu)=rc(pu),生成wc(pu);提升與校正:w(pS)=w(pS)+P(pS)wc(pu).

算法4 ParaRelaxBlockGS算法Step1 子進(jìn)程組GE求解AE(p)對應(yīng)的子系統(tǒng):若ddE=0,對AE(p)wE(p)=gE(p)采用BoomerAMG法作迭代;否則采用并行Jacobi法.Step2 子進(jìn)程組GE將wE(p)的數(shù)據(jù)分別發(fā)送給子進(jìn)程組Gα(α=R,I);子進(jìn)程組GR和GI分別接收wE(p)的數(shù)據(jù).Step3 子進(jìn)程組Gα(α=R,I)分別作以下步驟:(1)計算右端:fα(p)=gα-DαEwE;(2)若ddα=0,對Aα(p)wα(p)=fα(p)采用BoomerAMG法作迭代;否則采用并行Jacobi法.經(jīng)過以上三步,所得的wα(p),α=R,E,I即為所需的w(pS).

② 將并行數(shù)據(jù)I轉(zhuǎn)換為并行數(shù)據(jù)II.這一步的關(guān)鍵是生成并行矩陣A(pS)屬于G(k)的長方陣A(pS,k).由式(7)易知，只需以下拼接:(a) 在子進(jìn)程組GR上：將AR(p,i),VRE(p,i)拼接，可得A(pS,i)；(b) 在子進(jìn)程組GE上：將VER(p,i),AE(p,i),VEI(p,i)拼接，可得A(pS,i+m)；(c) 在子進(jìn)程組GI上：將VIE(p,i),AI(p,i)拼接，可得A(pS,i+2m).

假設(shè)上述數(shù)據(jù)轉(zhuǎn)換工作已經(jīng)完成，下面給出算法1的并行算法，它包含Setup和Solution兩個過程，首先給出Setup過程的算法(算法2).

接著介紹算法1的Solution過程(算法3)，為描述方便起見，對一個給定的向量v∈R3n，引入記號vα(p),α=R,E,I表示v(pS)限制在子進(jìn)程組Gα上的并行向量.

最后，給出算法3中ParaRelaxBlockGS的算法描述(算法4).

3　數(shù)值實驗

例1LARED-S程序中直角坐標(biāo)系下的三個典型三維三溫線性離散系統(tǒng)：S4-4、S5-5、S6-2，網(wǎng)格規(guī)模為125×125×125.收斂準(zhǔn)則取殘差向量范數(shù)下降7個量級.實驗環(huán)境：Mem 7 GB、CPU i5-6200U、2.3 GHz*4；編譯優(yōu)化參數(shù)：-O2.

如表1～2所示，對三維問題：(a) 在迭代次數(shù)上，自適應(yīng)PCTL法明顯少于Boomer AMG法，與PCTL法相近；(b) 在CPU時間上，相對Boomer AMG法和PCTL法，分別平均加速2.96和3.22倍.

例2LARED-S程序中球坐標(biāo)系下的6個典型二維三溫線性離散系統(tǒng)：M1-3、M32-1、M33-2(Mi-j：第i個時間步中的第j個非線性迭代)，網(wǎng)格規(guī)模為32 000×24.實驗環(huán)境：Mem 24 GB、CPU W5590、3.33 GHz*8；編譯優(yōu)化參數(shù)：-O2.

表1　三種PGMRES(30)法的迭代次數(shù)及CPU時間

表2　兩種PCTL的內(nèi)迭代次數(shù)

表3　兩種PGMRES(30)法的迭代次數(shù)

如表3所示：(a) 并行自適應(yīng)PCTL法對應(yīng)的迭代次數(shù)少于Boomer AMG法，特別對于M32-1和M32-3，前者不超過12次，而后者超過100次；(b) 迭代次數(shù)并不隨進(jìn)程數(shù)的增加而變動，表明并行自適應(yīng)PCTL法具有良好的算法可擴展性.

例3針對二維情形下典型的內(nèi)爆模型進(jìn)行數(shù)值模擬，考察前100個時間步(每個時間步固定5次非線性迭代)中涉及的線性代數(shù)系統(tǒng)的求解，其中網(wǎng)格規(guī)模為16 384×288.收斂準(zhǔn)則取殘差向量范數(shù)下降10個量級.實驗環(huán)境：某國產(chǎn)大規(guī)模并行計算機；編譯優(yōu)化參數(shù)：-O2.

如圖1所示，對于多核情形：(a) 并行自適應(yīng)PCTL法所需的迭代次數(shù)更少且更穩(wěn)定，與Boomer AMG法相比，CPU核數(shù)為384時加速1.75倍；(b) 并行自適應(yīng)PCTL法的迭代次數(shù)隨進(jìn)程數(shù)的增加基本不動，表明其良好的算法可擴展性.

[1]POMRANIIG G.The equations of radiation hydrodynamics[M]. Pergamon: Oxford, 1973.

[2]莫則堯, 張愛清, 曹小林，等. 多介質(zhì)輻射流體力學(xué)數(shù)值模擬中的并行計算研究[J].自然科學(xué)進(jìn)展,2006, 16(3): 287-292.

[3]聶存云, 舒適. 一類二維三溫輻射熱傳導(dǎo)方程組的對稱有限體格式[J].湘潭大學(xué)自然科學(xué)學(xué)報, 2004, 26(1): 17-22.

[4]YUE X,SHU S,XU X, et al. An adaptive combined preconditioner with applications in radiation diffusion equations[J].Commun Comput Phys, 2015,18(5): 1313-1335.

[5]YUE X,XU X,SHU S. JASMIN-based two-dimensional adaptive combined preconditioner for radiation diffusion equations in inertial fusion research[J].E Asian J Appl Math, 2017，7(3): 495-507.

[6]BALDWIN C,BROWN P N,FALGOUT R, et al. Iterative linear solvers in 2D radiation-hydrodynamics code: Methods and performance[J].J Comput Phys, 1999，154: 1-40.

[7]MO Z. Parallel adaptive solution for two dimensional 3-T energy equation on UG[J].Comput Visual Sci, 2006，9: 165-174.

[8]JIANG J,HUANG Y,SHU S, et al. Some new discretiztion and adaptation and multigrid methods for 2-D 3-T diffusion equations[J].J Comput Phys, 2007，224(1): 168-181.

[9]郭美珍, 江軍, 舒適. 一種二維三溫輻射熱傳導(dǎo)方程組SFVEM格式的并行預(yù)條件子[J].湘潭大學(xué)自然科學(xué)學(xué)報,2007，29(3): 42-45.

[9]徐小文, 莫則堯, 安恒斌. 求解二維三溫輻射擴散方程組的一種代數(shù)兩層迭代方法[J].計算物理, 2009，26(1): 1-8.

[10]周志陽, 徐小文, 舒適,等. 二維三溫輻射擴散方程組兩層預(yù)條件子的自適應(yīng)求解[J].計算物理, 2012，29(4): 62-70.

[11]MO Z,ZHANG A,CAO X, et al. JASMIN: a parallel software infrastructure for scientific computing[J].Front Comput Sci,2010，4(4): 480-488.

[12]JESSEE J,FIVELAND W,HOWELL L, et al. An adaptive mesh refinement algorithm for the radiative transport equation[J].J Comput Phys, 1998，139: 380-398.

[13]HENSON V,YANG U. BoomerAMG: a parallel algebraic multigrid solver and preconditioner[J].Applied Numerical Mathematics,2002，41(1): 155-177.