三溫輻射擴(kuò)散問(wèn)題的一種并行自適應(yīng)PCTL預(yù)條件子*

岳孝強(qiáng),　周志陽(yáng),　徐小文,　舒　適

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

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

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

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

1　模型問(wèn)題及有限體積格式

考慮三溫RDEs

(1)

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

(2)

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

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

(3)

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

算法1 自適應(yīng)PCTL法:w=Bαg對(duì)給定的兩組參數(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),則作以下三步:前磨光過(guò)程:分別求解電子、輻射和離子方程子系統(tǒng)AEwE=gE,ARwR=gR-DREwE,AIwI=gI-DIEwE.具體做法:若Aα滿足式(6),則采用Jacobi迭代求解;否則采用AMG法求解,取最大迭代次數(shù)、迭代控制精度分別為μα、ε.求解粗水平子系統(tǒng):(PTAP)wc=PT(g-Aw),其中P由式(3)計(jì)算.具體做法:采用AMG法求解,取最大迭代次數(shù)、迭代控制精度分別為μc、ε.粗水平校正:w=w+Pwc.

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

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

(7)

為下標(biāo)相匹配，以下用記號(hào)GR、GE、GI來(lái)分別表示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)程分組策略的實(shí)現(xiàn)，需要對(duì)數(shù)據(jù)式(8)進(jìn)行轉(zhuǎn)換，此時(shí)需指定分劃數(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)系.以下分兩步來(lái)完成上述并行數(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)的長(zhǎng)方陣.(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過(guò)程Step1 將進(jìn)程組G分成子進(jìn)程組GR、GE、GI.Step2 生成弱耦合條件和強(qiáng)對(duì)角占優(yōu)條件是否成立的標(biāo)志變量.具體為子進(jìn)程組Gα(α=R,I):利用Aα(p),按式(6),生成強(qiáng)對(duì)角占優(yōu)標(biāo)志標(biāo)量ddα;子進(jìn)程組GE:利用AE(p)和VER(p),VEI(p),按式(4)至式(6),生成強(qiáng)對(duì)角占優(yōu)標(biāo)志變量ddE和弱耦合條件標(biāo)志變量wcE,并將wcE廣播給其他兩個(gè)子進(jìn)程組.Step3 若wcE=1,則子進(jìn)程組Gα(α=R,E,I)分別作以下步驟:若ddα=0,則對(duì)Aα(p)調(diào)用BoomerAMG法[13]的Setup算法.Step4 若wcE=0,則執(zhí)行以下步驟4.1 三個(gè)子進(jìn)程組分別作以下步驟(1)子進(jìn)程組Gα(α=R,I):①對(duì)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,則對(duì)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 對(duì)Ac(pu)調(diào)用BoomerAMG法的Setup算法.

算法3 并行自適應(yīng)PCTL法的Solution過(guò)程:w=Bαg若wcE=1,則子進(jìn)程組Gα(α=R,I)分別作以下步驟,生成w(pS):以0為初值,對(duì)Aα(p)wα(p)=gα(p)作να次迭代,其中:當(dāng)ddα=1時(shí)采用并行Jacobi法作迭代,否則采用BoomerAMG法的Solve算法作迭代.若wcE=0,則執(zhí)行以下步驟:調(diào)用ParaRelaxBlockGS算法,作一次塊Gauss-Seidel磨光,生成w(pS);計(jì)算殘量: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)對(duì)應(yīng)的子系統(tǒng):若ddE=0,對(duì)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)計(jì)算右端:fα(p)=gα-DαEwE;(2)若ddα=0,對(duì)Aα(p)wα(p)=fα(p)采用BoomerAMG法作迭代;否則采用并行Jacobi法.經(jīng)過(guò)以上三步,所得的wα(p),α=R,E,I即為所需的w(pS).

② 將并行數(shù)據(jù)I轉(zhuǎn)換為并行數(shù)據(jù)II.這一步的關(guān)鍵是生成并行矩陣A(pS)屬于G(k)的長(zhǎng)方陣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兩個(gè)過(guò)程，首先給出Setup過(guò)程的算法(算法2).

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

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

3　數(shù)值實(shí)驗(yàn)

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

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

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

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

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

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

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

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

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

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

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

[3]聶存云, 舒適. 一類二維三溫輻射熱傳導(dǎo)方程組的對(duì)稱有限體格式[J].湘潭大學(xué)自然科學(xué)學(xué)報(bào), 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é)報(bào),2007，29(3): 42-45.

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

[10]周志陽(yáng), 徐小文, 舒適,等. 二維三溫輻射擴(kuò)散方程組兩層預(yù)條件子的自適應(yīng)求解[J].計(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.