球谐函数展开快速算法及其并行算法研究

球谐函数展开快速算法及其并行算法研究球谐函数展开快速算法及其并行算法研究国防学技术大学研究生院分类号 TP3 12 学号09060065,,, 密级公开工学硕士学位论文球谐函数展开快速算法及其并行算法研究硕士生姓名王翔学科专业计算机科学与技术研究方向计算机应用技术指导教师宋君强研究员国防科学技术大学研究生院二〇一一年十一月Research on Parallel Algorithms of the Fast Algorithm for Spherical Harmonic ExpansionsCandidate Wang XiangSupervisor Prof Song JunqiangA thesisSubmitted in partial fulfillment of the requirementsfor the degree of Master of Engineeringin Computer Science and TechnologyGraduate School of National University of Defense Technology ChangshaHunancomaNovember 2011国防科学技术大学研究生院工学硕士学位论文目录摘要 iABSTRACT ii第一章引言 111 数值天气预报与谱模式 1com 谱模式的发展与现状 1com 球谐函数与Silberman 方法2com 变换法4com 谱模式的优缺点 512 球谐函数 6com 球谐函数的推导 6com 截断问题 6com 球谐函数展开算法的发展 713 GPU 通用计算与数值天气预报 914 本文研究内容 1015 论文结构 11第二章 CPUGPU 异构集群研究 1221 GPU 发展的历史与现状 1322 CPUGPU 异构集群 14com CPU 与GPU 的性能对比 14com 典型CPUGPU 异构集群 1623 并行编程环境 19com MPI 并行程序设计原理 19com CUDA 架构与编程模型 20com MPICUDA 混合编程模型研究 22第三章球谐函数展开快速算法研究 2431 球谐函数展开简介 2432 Rokhlin-Tygert 球谐函数展开快速算法 25com SerSTEP1函数分解 25com SerSTEP2快速Fourier 变换 25com1 余弦变换25com2 正弦变换27第I 页国防科学技术大学研究生院工学硕士学位论文com SerSTEP3低阶连带Legendre 多项式展开系数 27 com SerSTEP4任意阶连带Legendre 多项式展开系数 29 33 0 阶Legendre 多项式展开系数的计算 31com 直接法 31com 变换法 32com 两种方法比较 3434 标准特征向量矩阵的计算 3535 RT 算法串行程序实现 36第四章 RT 算法的MPI 并行程序设计 3841 RT 算法的并行策略 3842 RT 并行算法描述 4043 RT 并行算法的MPI 实现 42com 程序结构 43com 通信机制 44com 同步策略 4544 RT 算法的MPI 并行程序实验及分析 46 com 实验方案I 47com 实验方案II 48com 实验方案III 49第五章 RT 算法的MPICUDA 并行程序设计 51 51 GPU 在RT 算法中的加速策略 51com GPU 在RT 串行算法中的加速策略 51 com GPU 在RT 并行算法中的加速策略 5152 库函数 52com CuFFT 52com LAPACK 53com CULA 55com CulaBLAS 5553 RT 算法的异构混合并行实现 56com RT 算法串行程序的CUDA 加速 56com RT 算法MPICUDA 异构并行实现 5754 RT 算法的异构并行程序实验 58com 实验方案I 59com 实验方案II 59第II 页国防科学技术大学研究生院工学硕士学位论文com 实验方案III 60com 实验方案IV 62com 实验方案V 62第六章结论与展望 64致谢 65参考文献 67作者在学期间取得的学术成果 71附录A 公式推导 72第III 页国防科学技术大学研究生院工学硕士学位论文表目录表31 四个系数矩阵 31表32 求0 阶展开系数两种方法比较 35表33 求标准特征向量矩阵四种方法比较 36表34 串行程序在不同规模下各部分时间对比 36表41 MPI 版RT 并行程序文件列表 43表42 MPI 版RT 并行程序通信函数列表 45表43 各进程在示例条件下所得行号 46表44 使用48 进程实验方案 47表45 使用24 进程实验方案 48表46 三种进程数实验对比 49表51 CuFFT 函数参数说明 53表52 CLAPACK 函数参数说明 54表53 culaSsyev 参数说明 55表54 culaSgemv 参数说明 56表55 CUDA 加速RT 串行程序文件列表 57表56 MPICUDA 异构并行RT 算法文件列表 58表57 SerSHE 与CUDASerSHE 实验对比 59表58 MPIParaSHE 在三种进程规模下的实验对比 60表59 MPICUDAParaSHE 程序64 节点实验对比 61表510 MPICUDAParaSHE 程序512 节点实验对比 62表511 MPIParaSHE 与MPICUDAParaSHE 程序对比实验 63 第 IV 页国防科学技术大学研究生院工学硕士学位论文图目录[1]图11 两种波数截断方式 7图12 北极附近的风速向量 8[27]图21 CPU 与GPU 的结构对比 12[27]图22 CPU 和GPU 峰值浮点计算能力比较 15[27]图23 CPU 和GPU 的带宽比较 16[36]图24 DEGIMA 的互连网络 17[55]图25 天河一号组织结构图 18图26 CUDA 的软件系统结构 21图27 CUDA 编程模型 22图28 MPICUDA 并行计算模型 23图31 Rokhlin-Tygert 球鞋函数展开快速算法 26图41 RT 并行算法流程 42图42 MPI 版RT 并行程序中文件间的关系 44图43 ParaSTEP4 程序示例 45图44 划分通信组程序示例 46图45 两种方案中的对比 48图46 求特征向量矩阵时间增长 49图51 三种进程规模的并行效率 60图52 MPICUDAParaSHE 在两种进程规模下加速比变化 63 第 V 页国防科学技术大学研究生院工学硕士学位论文第 VI 页国防科学技术大学研究生院工学硕士学位论文摘要计算机技术的迅速发展促进了数值天气预报模式的不断推进并使偏微分方程离散化相关的数值方法变得更加完善和精确从而减小了求偏微分方程数值解的误差谱变换方法采用球谐函数的有限截断展开来逼近球面某垂直层上的状态变量充分利用变量在物理空间和谱空间之间的适当变换有效降低偏微分方程中以非线性项为主的计算量进而获得高精度的计算解其中物理空间由状态变量在经纬度格点的直积组成谱空间由状态变量对应的截断球谐函数展开系数组成谱模式具有精度高稳定性好等优点但其运算量和存储量较大球谐函数展开作为一种数学方法被广泛应用于数值天气预报模式地球物理化学物理偏微分方程数值解等众多领域球谐函数是球面纬向微分算子以及球面Laplace 算子的特征函数基于球谐函数的谱变换是谱模式计算中的核心问题涉及 Fourier 变换和 Legendre 变换其中 Fourier 变换应用于经度方向Legendre 变换应用于纬度方向本文基于Rokhlin-Tygert 球谐函数展开快速算法深入研究了算法的并行化问题主要内容包括1分析数值天气预报谱模式以及球谐函数的特点研究了并行计算理论和技术包括MPI 并行CUDA 并行MPICUDA 异构并行2 深入研究并补充完善了Rokhlin-Tygert 球谐函数展开快速算法增加了0 阶标准Legendre 多项式展开系数的计算步骤3 针对球谐函数展开快速算法设计了RT 并行算法4 实现了RT 算法的MPI 并行程序并分别在某五万亿次集群和天河一号超级计算机上进行测试其结果均显示在三角截断波数M1023 时并行效率可达87 以上5 实现了RT 算法的MPICUDA 异构并行程序调用了CULA 库函数天河一号上的实验结果显示CUDA 的加速效果在M1023 时逐渐显著单节点双进程布局的并行效率随的增大而增长并超过单节点单进程在MM 4095时使用512 节点1024 进程可获得超过6700 的加速比主题词球谐函数展开并行算法GPU 通用计算CULAMPICUDA第 i 页国防科学技术大学研究生院工学硕士学位论文ABSTRACTThe continuous progress of Numerical Weather Pridiction NWP model is promoted by the rapic development of computer technology which makes moreconsummate and accurate the numerical method of discretization of continuous partialdifferential equations PDE and therefore decreases the error of numerical solutions ofPDE Spectral transform method means using spherical harmonic expansions SHEwith limited truncations to approximate some variable of a specific vertical layer ofsphere surface taking advantage of appropriate transformation of the variable betweenphysical and spectral space to gain accurate solutions and decrease the computationwhich mainly attributes to the non-linear terms of PDE The physical space is formed byCartesian products of the variable on the points of longitudinal-latitudinal grid whilethe spectral space is formed by the corresponding coefficients of truncated SHE of thevariable Spectral model is advanced with its high accuracy and stability while with thedrawback of large mass of computation and storageAs a mathematical method SHE is implemented in many domainssuch as NWPmodel Geophysics Chemical physics numerical solutions of PDE etc Sphericalharmonics are the eigenfunctions of latitudinal differential operators and Laplaceoperators on sphere surface Based on spherical harmonics thespectral transform is thekernel of the computation of spectral model consisting of Fourier Transform andLegendre Transform which is applied in longitudinal and latitudinal directionrespectivelyBased on Rokhlin-Tygert s fast algorithms for SHE RT Algorithm this dissertation deeply investigates the parallelization of the algorithm mainly including1 Analyze the characteristics of the spectral model of NWP and sphericalharmonics and study the theories and techniques of parallel computing including MPIparallelism CUDA parallelism MPICUDA heterogeneous parallelism2 Deeply study and consummate Rokhlin-Tygert s fast algorithms for SHE towhich adding a procedure for computing the coefficients of 0-order normalizedLegendre functions3 Design RT parallel algorithm according to the fast algorithms for SHE4 Implement MPI parallel program for RT Algorithm which was tested on the5-Terascale Clusters and Tianhe-1A supercomputer respectively the results of whichshow that when the wave number of triangular truncations M 1023 the parallelefficiency of the program would reach over 875 Implement MPICUDA heterogeneous parallel program for RT Algorithmusing CULA library functions The result of the test on Tianhe-1A shows that firstlyCUDA acceleration would be effective gradually when M 1023 secondly the第 ii 页国防科学技术大学研究生院工学硕士学位论文parallel efficiency of 2-process-per-node distribution would increase as M grows andeventually transcend that of 1-process-per-node distribution thirdly when M 4095the speedup of using 1024 processes on 512 nodes could reach over 6700Key Words Spherical Harmonic Expansions Parallel Algorithm GPUGeneral Computing CULA MPICUDA第 iii 页国防科学技术大学研究生院工学硕士学位论文第 iv 页国防科学技术大学研究生院工学硕士学位论文第一章引言11 数值天气预报与谱模式com 谱模式的发展与现状1922 年Richardson 首次提出的关于数值天气预报 NumericalWeatherPrediction 的设想在1950 年由Charney 等人在普林斯顿的高等研究所的ENIAC[1]计算机上实现他们应用正压过滤模式第一次成功地制作了24 小时的数值预报全球或半球范围的的大气数值模式的进展一方面表现在逐步减少模式方程组中所取的各种近似改进模式中所包含的各种物理过程的参数化方案使之能够更准确真实全面地反映实际大气中的各种动力和热力过程另一方面也表现在模式的水平和垂直分辨率在不断地提高数值天气预报模式的进展与计算机技术的迅速发展紧密地联系在一起而用于将连续偏微分方程离散化的数值方法也变得更加完善和精确从而减小求偏微分方程数值解的误差数值天气预报和大气环流数值模拟的实践已经证明有限差分法是求解流体动力学方程组的一种有效方法谱方法在数值天气预报模式中的应用可上溯至 1942 年当时前苏联学者Blinova[2]提出利用球谐函数展开Spherical Harmonic Expansions 的线性化模式来制作长期天气预报随后Silberman[2]在 1954 年借助于球谐函数求解了无辐散的正压涡度方程Lorenz[3]在 1960 年证明了正压无辐散流动的谱截断方程与其对应的精确微分方程一样具有总动能和总涡度平方守恒的性质Platzman[4]于同年指出这一性质可以自动地消除数值计算的非线性不稳定而这在当时是以有限差分法为基础的模式网格点模式中存在的最主要的困难之一1965 年 Cooley[5]等人提出的快速傅里叶变换Fast Fourier Transform FFT 使得谱方法的计算效率大大提高1966 comrt[6]说明了如何将展开函数分解为较简单的形式以解决计算机内存不足的问题1970 年 Eliasen[7]等人和 Orszag[8]分别独立提出变换法计算二次型的非线性项以取代相互作用系数法使得?谱模式所需的计算量与计算机内存容量大为降低?更容易将各种物理过程直接引入谱模式中变换法的提出为谱模式的发展和应用铺平了道路[9]我国对谱模式的研究始于20 世纪70 年代初郑庆林从 1974 年开始利用北半球三层准地转谱模式进行了较多个例的预报实验其72 小时内的预报效果是比[10]较好的他又于1976 年做出了北半球七层初始方程谱模式的个例预报该预报[11]模式考虑了地形和动量热量的水平扩散等物理过程1989 年郑庆林将该模第 1 页国防科学技术大学研究生院工学硕士学位论文式发展为全球七层大气环流谱模式并进行了30 天的长期数值天气预报试验该模式包含了较完整的物理过程方程组求解方案能有效地克服散度方程以及σ 坐标系中在大地形附近计算气压梯度力项时所存在的大量之间的小差问题其非线性项的谱计算方法有优越性且在整个30 天内模式的预报误差都比对应的持续性误[12]差小1991 年韩卫清纪立人通过一个单向嵌套的有限区域浅水波方程谱模式的具体方案预报试验验证了谱模式对于边界条件的敏感性并提出了一种有[13]助于谱模式中对陡峭地形处理的地形区计算方案1996 年陈嘉滨舒静君提出用热力学变量比容α代替温度变量T 并将地面气压本身作为展谱预报变量解决了模式大气总质量不守恒的问题且维持气压梯度项为二阶项不出现混淆误[14]差谷湘潜于1998 年推导出大气自忆方程的谱形式采用T42L9 谱模式为动力核建立了全球自忆T42 模式对500hPa 高度场的预报试验结果将月平均环流距平相关系数提高135月平均环流均方根误差降低294com 球谐函数与Silberman 方法在球坐标中相对涡度为1 vu cos11a cos无辐散流场可用流函数表示为1ua121va cos涡度表示为2 1 1 22 2 cos 13a cos cos则正压涡度方程可写成[56] 2 1 2 2 22 2 14t a cos a函数Ym Pm sin eim 15l l是方程2 2 m m 16a Y l l 1 Y 0l l的解其中P m sin 是连带Legendre 多项式Associated Legendre Functionsl第 2 页国防科学技术大学研究生院工学硕士学位论文进而流函数可按球谐函数展开为J m J2 m m17t a t Yl lmJ l m从而涡度按球谐函数展开为J m J2 m m18l l 1 t Yl lmJ l m将公式17 和18 代入14交换非线性平流项的求和附标 m l 和 m l 1 1 2 2二次写出该项的表达式两式相加再除以2 可得J m J dm J m J il m m ml l 1 Y 2i m Y l ll l l 2 2mJ l m dt m J l m 2cos m1 l1 m2 l2 19m mdP 2 dP 1l l i m mm1 m2 m1 2 m2 1 1 21 l l 1 m P m P e1 1 l1 l2 1 l1 2 l2d d2 2公式19 各项乘以Ym cos 并积分 d d 可得l2 0dm 2im il m m1 m2L l l l m m m 110l l1 l2 1 2 1 2dt l l 1 2m1 l1 m2 l2其中相互作用系数l l 1 l l 1 2 dPm21 l2 2 1 1 P m m P m1 22 l l 1 2 l 1 l1 ddPm1111L l l l m m m m2 l11 2 1 2 m P d m m m2 l 1 22 d0 m m m1 2m1 m2 m表示和两个分量相互作用对于所产生的影响l1 l2 l由公式111可见计算一个相互作用系数所需计算复杂性为O J 下面计算相互作用系数的个数l l 2对于和两层求和由于二者求和次数均为则这两层总共为J J1 2m m m m对于和两层求和由于固定且仅当二者之和为时相互作用系1 2数不为0因此需要计算的相互作用系数的个数为O mm 的取值范围为 J J 则个数为O J可得相互作用系数的个数为J 2 O m O J O J 4因此使用Silberman 方法计算所有的相互作用系数其计算复杂性为O J 5 另一第 3 页国防科学技术大学研究生院工学硕士学位论文方面要存储所有的相互作用系数包括0 也要求O J 5 的存储量所以计算2机的运算速度和存储器大小直接影响了使用 Silberman 方法求解预报问题的精度com 变换法1970 年Eliasen 等人和Orszag 分别独立提出变换法用于计算非线性平[56]流项大大降低了谱方法的计算复杂性令1 2 2J 1122a cos在计算时不直接使用相互作用系数而是用和的展开系数计算网格点J上的涡度平流J 再由J 求出其球谐函数展开即非线性平流项是经过由谱到网格再由网格到谱的两次变换所得变换法中与公式110对应的谱方程为dm 2im 1l m mJ113l ldt l l 1 l l 1设为格点坐标公式112中各项分别定义为j k1 2 J m J 1 m YmF a l1 i k lcos j k mJ l m cosj j kJ m J2 m mF iml l 1 Y2 i k l l i kj k mJ l m114J m J2 m mF a im Y3 i k l l i kj k mJ l mJ m J l l 1 m1 2 m YF l4 i k lcos j k mJ l m cosj j k1 Ym其中 l 可由递推公式cos2 Ym Ym m m m mcos l cos l l Y l 1 Y 115l1 l1 ll 1sin计算其中1 2l2 m2m116l 4l2 1将公式114代入112可得第 4 页国防科学技术大学研究生院工学硕士学位论文1J F F F F 117i k 2 1 i k 2 i k 3 i k 4 i kaJ m由公式117计算格点值J 再进行球谐函数展开可得最后对方i k l程113进行时间积分由公式114 117可见计算一个格点值所需计算复杂性为O J 2 所以对于球面上的一维函数格点数为O J 从而计算复杂性由Silberman 方法的O J 5 降为O J 3 对于二维球面网格格点数为O J 2 因此计算一个固定阶数J m O J 2的非线性平流项所需存储量为则对于O J 个非线性平流项的计算所l需存储量为O J 3 远少于Silberman 方法的O J 5com 谱模式的优缺点谱变换方法的基本做法是用球谐函数的有限截断展开逼近球面某垂直层上的状态变量充分利用物理空间和谱空间之间的适当转换降低偏微分方程中以非线性项为主计算量获得高精度计算解其中物理空间由状态变量在经纬度格点的直积组成谱空间由状态变量对应的截断球谐函数展开系数组成与网格点模[1]式相比谱模式的优点主要表现在以下六个方面1计算空间微商精确使谱方法估计的位相速度比一般差分法更准确2 对于二次型的非线性项计算消除了非线性混淆现象可避免由此引起的计算不稳定3 用谱方法求解球坐标下的控制方程组不需要像有限差分法那样对球面网格中的极点做特殊处理因而特别适合于全球或半球模式4 在谱模式中易于应用半隐式时间积分方案其计算程序比网格点模式简单5 谱方法能自动并彻底滤去短波效果比一般差分法中用平滑算子好6 由于在全球模式中通常选择球谐函数作为谱展开式的基函数而球谐函2p数正好是球面上的Laplace 算子的特征函数所以在模式中计算 p 为正整数型的水平扩散项非常方便基于同样的理由用谱方法解Poisson 方程或Helmholtz方程也特别方便不需要进行迭代[1]另外谱方法的缺点体现在以下四个方面1运算量和存储量较大当模式的水平分辨率提高时谱方法的计算量比格点法增加得更快2 对分布不太连续的物理量容易发生Gibbs 现象需要用较多的谱分量第 5 页国防科学技术大学研究生院工学硕士学位论文才能表示3 当m?0 时在高纬度连带Legendre 函数P的值很小用球谐函mn数展开地形高度误差较大4 制作有限区域预报或套网格预报不如差分法灵活方便12 球谐函数com 球谐函数的推导在球坐标系中Laplace 方程为21 u 1 u 1 u22 r 2 cos 2 2 2 0 118r r r r cos r cos其中为为到球心的距离为纬度为经度用分离变量法求解方程15时r可令u r R r Y 119带入方程 118可得2r R 2rR R 021202 1 Y 1 YY Y cos Y 02 2cos cos[1]其中第2 个公式称为球面调和方程满足球面调和方程的有连续二阶导数的有界解称为球谐函数如公式15com 截断问题在公式1718114中所有的内层求和上下界之差均为即对Jm于每一个阶数球谐函数展开次数相同其中为纬向波数l m 为经向波m[56]数在谱模式中应用最普遍的两种波数截断方式称为三角形截断和平行四边形截断如图11 所示其中三角形截断为J JmJ l m所有展开的波都限制在波数平面上的一个三角形区域内平行四边形截断为J m JmJ l m[1]展开的波都限制在波数平面上的一个平行四边形区域内第 6 页国防科学技术大学研究生院工学硕士学位论文l lmma 三角形截断b 平行四边形截断图11 两种波数截断方式[1]早期的谱模式几乎都使用平行四边形截断com 节中的Silberman com节中的变换法但是大多数后来投入正式业务使用的谱模式使用三角形截断如法国日本欧洲中期天气预报中心的模式选择截断方式的原则应最充分利用给定的自由度来获得最精确的计算结果或在一定的精度要求下花费最小的计算代价在早期的谱模式中应用平行四边形截断有其历史原因一方面平行四边形截断的谱分量可以存储在一个二维数组中程序编写较容易另一方面Ellsaesser[57]证明了对于甚低分辨率的模式平行四边形截断可以使500hPa 旋转流场的动能保留最大的方差但是1972年Baer[58]在分析了1969年1-2 月北半球700500300250 200hPa 五层的动能谱厚倾向于在中等分辨率的模式中应用三角形截断三角形截断的一个突出优点是具有各向同性所得到的分辨率在整个球面上是均匀的即无论将球坐标系的极轴定在何处截断后的有限项展开式所表示的场都不会变化而平行四边形截断的有限项展开式所表示的场仅仅对于围绕地球极轴的坐标旋转变换保持不变即分辨率仅在东西向是均匀的在南北向既不一定与东西向分辨率相等还随地理位置变化从计算的角度看三角形截断要比同样自由度的平行四边形截断的计算量稍微节省首先这是因为平行四边形截断包含了较短的波分量所用时间步长就要相应地取短其次其所采用的沿经圈方向的计算格点数要比三角形截断多 20左右这两个原因似的平行四边形截断的计算代价要比同样自由度的三角形截断多25左右com 球谐函数展开算法的发展对于全球或半球的大气模式宜选用球面曲线坐标系而对于解球面坐标系[1]的控制方程组若采用有限差分法则存在极区问题如图12 所示在北极V V点附近从点A 趋近点P 时风速极限值为北风从点C 趋近点P 时为南风P。