岳孝强,　周志阳,　徐小文,　舒　适

(1.湘潭大学 数学与计算科学学院,湖南 湘潭 411105;2.北京应用物理与计算数学研究所,北京 100094)

在惯性约束聚变(ICF)等高能量密度物理领域的数值模拟中，辐射流体力学方程组的求解是非常重要的环节.该方程组由描述流体运动与辐射传输的两类方程耦合而成，其中后者在扩散近似下由三温辐射扩散方程组(RDEs)刻画[1].由于物理量间的复杂耦合等特点，三温RDEs离散系统的条件数极差，致使求解时间通常占整体模拟时间的70%以上[2-3].因此，研究三温RDEs的高效并行解法器是一项具有重要理论意义与实用价值且难度很大的课题.预条件Krylov子空间迭代法是求解三温RDEs离散系统的主要方法，其中影响算法效率的关键是预条件子的构造，目前已有很多研究工作[4-10].值得一提的是，[9][10]针对Euler网格构造了基于物理量粗化及耦合强度的PCTL法及其自适应方法(但限于串行实现).这些预条件算法(主要是ILU类、多重网格类以及两者的有机组合)和理论主要是针对二维结构和协调非结构网格，仍需不断完善与发展，特别是关于算法的并行设计以及可扩展性仍需进一步提升.

自适应结构网格应用支撑框架(JASMIN)以其数据结构、高效算法以及并行通信等高度封装的特性而成为大规模实际应用需求的一个重要研发平台[11].基于细化模式的隐式时间积分算法是求解ICF内爆过程中多物理耦合问题的常用模式[12].在该模式下，求解网格片层次结构上离散系统的主要工作量体现在每个网格层上离散系统的求解.因此，如何在不含悬点的分片结构网格上，为三温RDEs的大规模离散系统设计出基于自适应PCTL预条件子的高效并行PGMRES解法器，将为解决基于细化模式的ICF瓶颈问题起到积极作用.

本文针对三温辐射扩散问题的全隐有限体积格式，将自适应PCTL法推广到三维情形，给出已有的二维自适应PCTL法在JASMIN平台下的并行实现算法(涉及进程分组及数据转换)，并通过二维和三维典型算例验证了新并行解法器的高效性与良好的算法可扩展性.

1　模型问题及有限体积格式

考虑三温RDEs

(1)

(1)是非线性偏微分方程组，利用向后Euler方法进行时间离散、“凝固系数法”线性化，再采用全隐有限体积格式，可得大规模离散系统Au=f，其中A∈R3n×3n，n为网格规模.假定自由度按照光子、电子、离子进行排列，则A可写为块三对角结构

(2)

2　一种基于JASMIN平台的并行自适应PCTL法

2.1　一种自适应PCTL法

(3)

需要注意的是，自适应PCTL法的粗化过程是基于物理量进行的，算法描述.

算法1 自适应PCTL法:w=Bαg对给定的两组参数θ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)则调用να次Jacobi迭代求解;否则采用AMG法求解.若不满足式(4)或式(5),则作以下三步:前磨光过程:分别求解电子、辐射和离子方程子系统AEwE=gE,ARwR=gR-DREwE,AIwI=gI-DIEwE.具体做法:若Aα满足式(6),则采用Jacobi迭代求解;否则采用AMG法求解,取最大迭代次数、迭代控制精度分别为μα、ε.求解粗水平子系统:(PTAP)wc=PT(g-Aw),其中P由式(3)计算.具体做法:采用AMG法求解,取最大迭代次数、迭代控制精度分别为μc、ε.粗水平校正:w=w+Pwc.

2.2　自适应PCTL法的并行实现

由算法1可知：自适应PCTL预条件子涉及若干子系统的求解，且其中有些子系统可以独立求解.基于这个发现，本节将利用进程分组策略给出算法1的并行实现.为此，首先给出基于进程分组的分划数组的概念.

(7)

为下标相匹配，以下用记号GR、GE、GI来分别表示G1、G2、G3.

在JASMIN中，式(2)的并行数据为各子块的并行矩阵和交叉块的并行向量：

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

(8)

为便于算法1基于进程分组策略的实现，需要对数据式(8)进行转换，此时需指定分划数组p，这里我们将其取为p(j)=pu(3j),j=0,1,…,m.数据转换的目标是生成

(1) 并行数据I：各子进程组上的并行数据

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

式中：Aα(p)(α=R,E,I)表示子矩阵Aα在子进程组Gα上基于分划数组p的并行矩阵；Vαβ(p)(α,β=R,E,I)表示Vαβ在子进程组Gα上基于分划数组p的并行向量.

(2) 并行数据II：进程组G上基于分划数组pS的并行矩阵A(pS).

上述并行数据转换的本质是将矩阵和向量数据在各进程上进行分配，其关键是如何建立接口并行数据和目标并行数据之间的映射关系.以下分两步来完成上述并行数据转换.

① 将并行数据式(8)转换为并行数据I.(a) 将并行矩阵Aα(Pu)转为Aα(p),α=R,E,I.具体的数据分配方式为：将Aα(pu,3i)，Aα(pu,3i+1)，Aα(pu,3i+2)分配给进程Gα(i)，这里Aα(pu,k)表示并行矩阵Aα(pu)属于G(k)的长方阵.(b) 将并行向量Vαβ(pu)转为Vαβ(p)，α,β=R,E,I.具体的数据分配方式为：将Vαβ(pu,3i)，Vαβ(pu,3i+1)，Vαβ(pu,3i+2)分配给进程Gα(i)，这里Vαβ(pu,k)表示并行向量Vαβ(pu)属于G(k)的子向量.

算法2 并行自适应PCTL法的Setup过程Step1 将进程组G分成子进程组GR、GE、GI.Step2 生成弱耦合条件和强对角占优条件是否成立的标志变量.具体为子进程组Gα(α=R,I):利用Aα(p),按式(6),生成强对角占优标志标量ddα;子进程组GE:利用AE(p)和VER(p),VEI(p),按式(4)至式(6),生成强对角占优标志变量ddE和弱耦合条件标志变量wcE,并将wcE广播给其他两个子进程组.Step3 若wcE=1,则子进程组Gα(α=R,E,I)分别作以下步骤:若ddα=0,则对Aα(p)调用BoomerAMG法[13]的Setup算法.Step4 若wcE=0,则执行以下步骤4.1 三个子进程组分别作以下步骤(1)子进程组Gα(α=R,I):①对Aα(p)调用BoomerAMG法的Setup算法;②调用BoomerAMG法的Solve算法求解:Aα(p)pα(p)=-VαE(p),得pα(p).(2)子进程组GE:若ddE=0,则对AE(p)调用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)调用BoomerAMG法的Setup算法.

算法3 并行自适应PCTL法的Solution过程:w=Bαg若wcE=1,则子进程组Gα(α=R,I)分别作以下步骤,生成w(pS):以0为初值,对Aα(p)wα(p)=gα(p)作να次迭代,其中:当ddα=1时采用并行Jacobi法作迭代,否则采用BoomerAMG法的Solve算法作迭代.若wcE=0,则执行以下步骤:调用ParaRelaxBlockGS算法,作一次块Gauss-Seidel磨光,生成w(pS);计算残量:r(pS)=g(pS)-A(pS)w(pS);残量限制:rc(pu)=P(pS)Tr(pS);调用BoomerAMG法求解Ac(pu)wc(pu)=rc(pu),生成wc(pu);提升与校正:w(pS)=w(pS)+P(pS)wc(pu).

算法4 ParaRelaxBlockGS算法Step1 子进程组GE求解AE(p)对应的子系统:若ddE=0,对AE(p)wE(p)=gE(p)采用BoomerAMG法作迭代;否则采用并行Jacobi法.Step2 子进程组GE将wE(p)的数据分别发送给子进程组Gα(α=R,I);子进程组GR和GI分别接收wE(p)的数据.Step3 子进程组Gα(α=R,I)分别作以下步骤:(1)计算右端:fα(p)=gα-DαEwE;(2)若ddα=0,对Aα(p)wα(p)=fα(p)采用BoomerAMG法作迭代;否则采用并行Jacobi法.经过以上三步,所得的wα(p),α=R,E,I即为所需的w(pS).

② 将并行数据I转换为并行数据II.这一步的关键是生成并行矩阵A(pS)属于G(k)的长方阵A(pS,k).由式(7)易知，只需以下拼接:(a) 在子进程组GR上：将AR(p,i),VRE(p,i)拼接，可得A(pS,i)；(b) 在子进程组GE上：将VER(p,i),AE(p,i),VEI(p,i)拼接，可得A(pS,i+m)；(c) 在子进程组GI上：将VIE(p,i),AI(p,i)拼接，可得A(pS,i+2m).

假设上述数据转换工作已经完成，下面给出算法1的并行算法，它包含Setup和Solution两个过程，首先给出Setup过程的算法(算法2).

接着介绍算法1的Solution过程(算法3)，为描述方便起见，对一个给定的向量v∈R3n，引入记号vα(p),α=R,E,I表示v(pS)限制在子进程组Gα上的并行向量.

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

3　数值实验

例1LARED-S程序中直角坐标系下的三个典型三维三温线性离散系统：S4-4、S5-5、S6-2，网格规模为125×125×125.收敛准则取残差向量范数下降7个量级.实验环境：Mem 7 GB、CPU i5-6200U、2.3 GHz*4；编译优化参数：-O2.

如表1～2所示，对三维问题：(a) 在迭代次数上，自适应PCTL法明显少于Boomer AMG法，与PCTL法相近；(b) 在CPU时间上，相对Boomer AMG法和PCTL法，分别平均加速2.96和3.22倍.

例2LARED-S程序中球坐标系下的6个典型二维三温线性离散系统：M1-3、M32-1、M33-2(Mi-j：第i个时间步中的第j个非线性迭代)，网格规模为32 000×24.实验环境：Mem 24 GB、CPU W5590、3.33 GHz*8；编译优化参数：-O2.

表1　三种PGMRES(30)法的迭代次数及CPU时间

表2　两种PCTL的内迭代次数

表3　两种PGMRES(30)法的迭代次数

如表3所示：(a) 并行自适应PCTL法对应的迭代次数少于Boomer AMG法，特别对于M32-1和M32-3，前者不超过12次，而后者超过100次；(b) 迭代次数并不随进程数的增加而变动，表明并行自适应PCTL法具有良好的算法可扩展性.

例3针对二维情形下典型的内爆模型进行数值模拟，考察前100个时间步(每个时间步固定5次非线性迭代)中涉及的线性代数系统的求解，其中网格规模为16 384×288.收敛准则取残差向量范数下降10个量级.实验环境：某国产大规模并行计算机；编译优化参数：-O2.

如图1所示，对于多核情形：(a) 并行自适应PCTL法所需的迭代次数更少且更稳定，与Boomer AMG法相比，CPU核数为384时加速1.75倍；(b) 并行自适应PCTL法的迭代次数随进程数的增加基本不动，表明其良好的算法可扩展性.

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

[2]莫则尧, 张爱清, 曹小林，等. 多介质辐射流体力学数值模拟中的并行计算研究[J].自然科学进展,2006, 16(3): 287-292.

[3]聂存云, 舒适. 一类二维三温辐射热传导方程组的对称有限体格式[J].湘潭大学自然科学学报, 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]郭美珍, 江军, 舒适. 一种二维三温辐射热传导方程组SFVEM格式的并行预条件子[J].湘潭大学自然科学学报,2007，29(3): 42-45.

[9]徐小文, 莫则尧, 安恒斌. 求解二维三温辐射扩散方程组的一种代数两层迭代方法[J].计算物理, 2009，26(1): 1-8.

[10]周志阳, 徐小文, 舒适,等. 二维三温辐射扩散方程组两层预条件子的自适应求解[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.