引言

数值预报预测是提高气象预报预测能力的核心技术方法，数值模式是开展数值预报预测不可或缺的工具。高性能计算技术的应用在气象数值模式发展的历史上占有极其重要的地位，是气象数值预报预测业务发展的技术基础。几乎在任一时期，数值模式都使用了当时最先进的高性能计算机。数值模式向着更高分辨率、更复杂物理过程、集合预报、多模式耦合的方向发展，而模式分辨率越高，物理过程越精细，模式积分计算量越大，对计算能力的要求也越来越高^[1-3]。处理器的单核性能增长放缓，高性能计算系统性能的提升主要来自更大规模处理器数量，2010年前后至今，GPU(graphic processing units)、Intel Xeon Phi等众核处理器逐渐成为高性能计算系统的重要组成部件，高性能计算系统体系架构由多核CPU同构系统向异构系统发展^[4]。我国于2016年推出由全国产SW26010异构众核芯片构建的世界上第1台计算速度超过10亿亿次的高性能计算机神威·太湖之光系统^[5]，标志着我国高性能计算机的研制能力达到了国际领先水平。

体系架构的改变对各领域高性能计算软件带来明显影响，对于气象数值模式而言，由于具有成千上万行基于多核通用CPU开发的代码^[6-8]，向新型计算架构系统的迁移过渡尤为困难。为了使气象数值模式充分适应异构高性能计算体系架构，世界各大气象超算中心均持续开展相关研究。美国的WRF，GEOS-5，NIM/FV3，MPAS和欧洲的COSMO，IFS及日本的ASUCA等模式均基于NVIDIA GPU和Intel Xeon Phi系统开展了优化^[9-15]。美国CAM4-Chem模式在NVIDIA GPU进行移植，获得约3.8倍加速^[16]，我国BCCAGCM(BCC Atmospheric General Circulation Model)模式辐射过程在Intel KNC众核平台进行优化，获得约2.3倍性能提升^[17]。我国清华大学、中国科学院软件所、国家海洋局第一海洋研究所等单位也基于神威·太湖之光异构众核系统开展了数值模拟程序的性能优化^[18-19]。

中国气象局现有数值模式也存在这个问题，为了促进我国主要数值模式对新型计算架构的适应性应用，在神威·太湖之光系统开展BCCAGCM模式异构众核优化工作，为未来更高分辨率数值模式新型计算架构细粒度并行工作提供技术探索和积累。

1 神威·太湖之光超级计算系统简介 1.1 硬件系统

神威·太湖之光超级计算系统全部采用国产申威异构众核处理器SW26010构建，峰值运算性能125.44 PFLOPS，持续运算速度93.02 PFLOPS，性能功耗比6051.13 MFLOPS/W，内存总容量1024 TB，磁盘总容量20 PB。2016年开始，4次蝉联世界TOP500超级计算机榜单榜首^[20]。

SW26010异构众核处理器采用片上计算阵列集群和分布式共享存储相结合的异构众核体系结构。每个SW26010处理器集成了4个运算核组(core group, CG)，每个运算核组包含1个运算控制核心(称为主核，management processing element，MPE)和64个运算核心(称为从核，computing processing element，CPE)，64个从核构成1个从核阵列(CPE cluster)^[21]，架构如图 1所示。4个核组的物理空间统一编址，运算控制核心和运算核心均可访问芯片上的主存空间。主核是具有完整功能的64位RISC芯片，支持中断、内存管理、超标量和乱序执行等功能，可以进行通信、输入/输出、计算等操作。从核是具有有限功能的64位RISC芯片，提供聚合计算能力。每个从核有独立的、容量为64 KB的高速缓存(local data memory, LDM)，支持直接内存存取(direct memory access，DMA)的方式在主存和LDM之间传输数据。每个SW26010处理器峰值性能达到3.168 TFLOPS。

1.2 软件系统

神威·太湖之光系统采用定制的64位Linux操作系统，系统软件包括异构平台作业调度、资源管理、并行文件系统等系统软件。并行程序开发环境提供支持神威·太湖之光异构计算的编译、性能分析、并行调试等功能，支持主流编程语言及编程标准，包括C/C++, Fortran, MPI, OpenMP等。核组之间支持消息传递模型、核组内部支持共享内存模型，同时也支持加速线程等编程模型^[22-24]，软件环境如图 2所示。

1.3 神威·太湖之光OpenACC编程模型

OpenACC是2011年由美国Cray，PGI和NVI-DIA等公司联合推出的编程标准，可基于GPU等处理器加速科学计算代码，实现较为简单，具有可移植性，目前在异构编程中得到较为广泛的应用。OpenACC在串行C/C++或Fortran代码上添加编译制导语句，编译器将代码编译链接成并行程序，实现更高的性能。主要包括计算并行化、数据管理和运行环境等部分。由于SW26010芯片与GPU等异构计算平台在硬件体系结构存在一定差异，OpenACC在神威·太湖之光系统的实现和应用也有所不同，对标准OpenACC语义和功能进行了适当外延。神威·太湖之光OpenACC工具具有并行任务管理、异构代码管理和数据传输控制等功能，支持SWACC/SWAFORT编译器，编译过程如图 3所示。

1.3.1 执行模型

程序首先在主核启动运行，执行到由神威·太湖之光OpenACC编译制导语句标示的计算密集核心代码A，则将核心代码A和所需数据作为加速任务加载到从核阵列执行，执行完成后，将结果返回到主核，程序继续在主核运行，执行到下一段计算密集核心代码B，再将其加载到从核阵列执行，依次类推，将程序中所有计算密集核心代码加载到从核阵列执行，实现计算加速。执行过程如图 4所示。

1.3.2 存储模型

标准OpenACC支持主机内存和加速器内存分离的存储模型，如多数GPU平台。与此不同的是，SW26010芯片的主核和从核共享内存，即从核可直接访问主核数据空间，可称为共享内存的异构架构。由于直接访问主机内存导致性能损失，需要充分利用从核阵列中的64 KB LDM，提升数据访问效率。这也是神威·太湖之光OpenACC对OpenACC标准进行的主要功能扩展。在SW26010芯片中，数据空间主要包括主线程数据空间、加速线程私有空间和加速线程本地空间3类。存储模型见图 5。

由OpenACC执行模型和存储模型可知，计算加速效果收益主要来自于将可细粒度并行的计算密集代码加载到从核阵列并发执行。执行过程中，首先在从核阵列为所需加载数据分配空间，将代码和数据传输到从核阵列进行计算，再将计算结果传回主存，释放从核阵列数据空间。数据传输、控制管理等操作会增加额外负载开销，只有当从核阵列并行带来的收益大于负载开销，才可能获得一定性能提升。加速效果的获得与数据管理策略和并行算法设计密切相关，同时也与主核、从核自身能够提供的计算能力相关。

2 BCCAGCM模式移植 2.1 计算框架

BCCAGCM模式是中国气象局国家气候中心自主研发的全球大气环流模式，也是北京气候中心气候系统模式BCCCSM(Beijing Climate Center Climate System Model)的大气分量模式。该模式是主要短期气候预测业务模式^[25]，并成功应用于第5次国际气候变化影响评估报告(IPCC AR5)相关的气候模拟试验^[26]。模式主体采用Fortran90语言，部分采用Fortran77和C语言实现，包含近200个程序，约60万行代码。

BCCAGCM模式具有高度模块化结构如图 6所示。代码设计支持分布式和共享内存并行架构，包括：动力学求解大气动力方程及物理学近似求解网格中物理现象(如凝结过程、云、长短波辐射及湍流混合等)。在每个模拟周期内，动力学和物理学两个阶段交替进行。动力计算和物理计算在数据结构和并行策略方面均不相同，通过内部耦合将数据在两个阶段的数据结构之间进行相互转换。两个阶段分别针对自己的数据结构进行维护和优化^[27]。

2.2 移植策略

由于BCCAGCM模式基于通用CPU开发，向神威·太湖之光异构众核系统移植过程中，首先需要完成模式到主核的移植工作。通过修改编译器、编译选项，安装相应软件，对部分代码进行适应性修改。BCCAGCM模式可在申威主核并行执行，形成申威主核版本。为了获得更优的计算性能，需要对模式核心计算代码进行重构，使其可利用从核阵列进行加速，将每个核心段由主核运行方式扩展到主核与从核异构并行方式，并尽可能控制浮点计算误差。

异构并行的实现主要考虑两个问题：一方面要调整并行计算策略，另一方面要进行数据存储策略的适应性改写。BCCAGCM模式在通用CPU系统采用MPI与OpenMP混合并行策略，MPI用于节点间进程并行，OpenMP用于节点内线程并行。针对申威芯片主核和从核阵列的硬件架构，将原有的MPI与OpenMP并行方式改为MPI与OpenACC并行方式。调整部分代码的计算顺序和循环结构，获得足够的并行计算量，充分利用从核阵列的并行度，发挥其计算优势。由于每个从核的LDM具有较高的性能，由用户设计数据存取策略，控制数据在主存和LDM之间的传输，将适合的数据加载到LDM，提高已加载数据的可用性，降低数据的传输开销。

3 异构众核优化 3.1 性能分析

BCCAGCM模式具有计算构造复杂，代码量大的特点。首先对申威主核版本进行计算性能分析，选定计算耗时较多，程序算法可进行众核计算改造的核心段。设计程序测试算例，BCCAGCM模式选用T266(约50 km)水平分辨率，垂直方向共26层，全球水平网格为800×400，时间步长为180 s。在申威主核采用100 MPI任务，连续积分3 d，统计耗时排名，不考虑系统开销函数和标准数学函数，排名前几位的核心段及其占用总计算时间的百分比如表 1所示。

3.2 基本方法

采用从核阵列对程序进行计算加速，对程序代码进行修改，将核心代码由在主核执行加载到从核阵列执行，主要包含3个步骤：①在程序中寻找并行度高的核心循环段，核心循环中的代码须为串行执行，不能包含消息通信功能的代码。②在核心循环代码前面增加编译制导语句。确定核心代码中需要私有化的变量，并添加私有化子句。③编译运行，确定运行结果是否正确。在从核阵列执行的代码要获得更好的计算性能，需要对其进一步优化，主要考虑几个方面：充分利用从核LDM空间，关键数据应尽量存放在LDM，尽量缩小数据传输和计算的开销比例，合理设置循环分块大小，提高数据传输效率。

3.3 并行优化方案

由于大气模式物理过程采用单柱计算，每个柱内计算为串行，各柱体之间没有数据依赖，计算区域各柱体计算可并行执行。对于计算逻辑较为简单的循环可直接添加制导语句，将循环内各柱体计算部分分散到从核进行并行计算，达到计算加速的目的。对于涉及到中间变量数组数据量较大的情况，考虑从核本地存储空间受限，需反复进行多次数据拷贝传输，产生的数据通信开销影响从核加速效果，有时甚至出现负的性能加速。针对这种情况，可对外围较大函数进行分析，通过子函数进行循环内移的方式达到子函数分块，计算独立，再添加制导语句实现从核并行加速的目的。针对神威·太湖之光芯片主从核架构特点，充分利用LDM空间，缓解数据局存压力，降低通信开销，对于经过循环内移后的部分子函数循环方式进行调整，扩大循环计算量，提高并行度，使从核承担的计算任务与本地存储相适应。将外围循环变量提至循环外，并对原代码进行修改，必要时对数组进行升维。

3.4 核心段优化

对性能分析结果计算占比较高的重力波过程、勒让德逆变换系数计算和大尺度凝结降水过程等核心段分别进行众核加速优化，选择具有代表性的两部分进行介绍。

重力波过程计算耗时最多，对代码进行分析可知，该函数接口有较多的参数数组，传递内部操作所需的输入输出数据，函数内有部分临时数组，较多嵌套的循环迭代，大部分存在并行性，可用从核进行加速，但部分内部数组存在离散访存现象。首先，对传入数组和临时数组进行转置处理，将离散访存改为顺序访存。然后对函数内每段循环嵌套迭代进行计时，按运行时间长短依次采用从核进行优化。对部分函数内部数组，通过设置临时数组，修改声明与使用方法，消除离散访存，优化LDM使用。使用编译制导语句将较小循环合并，并行调度运行。对函数中多段循环进行优化后，重力波过程整体计算效率为未优化计算的8倍左右。

对大尺度凝结降水过程代码进行分析，结构比较复杂，循环内包括多个循环分块和多个函数调用，外层循环并行度较低，不能充分发挥从核并行计算能力。考虑物理过程计算模型为单柱计算，针对每个柱体进行多个物理过程计算，将物理过程进行循环内移。一方面有利于扩大从核任务规模，另一方面有利于对复杂的物理过程函数进行拆分。循环内移之后，添加制导语句，并采用分块操作，提高从核并行效率。进一步对从核访存性能进行优化，将全局变量的使用通过传递函数参数的方式来实现，避免位于程序代码最内层循环数据重复复制的现象。经过优化后，大尺度凝结降水过程整体计算效率为未优化计算的9倍。

3.5 优化效果

经过对部分动力及物理计算密集型核心段代码分析，采用循环内移、主从核访存优化、数据重构等方法进行优化，尽可能减少访存开销，优化核组内多从核执行加速性能。确保计算结果合理的前提下，各核心段分别获得一定性能提升，基本达到未优化计算的3倍左右，最高可达约14倍。各核心段函数加速及时间占比如图 7所示，计算量占比最高的重力波过程经优化后效率提升为原来的8倍，动力部分占比较高的勒让德逆变换系数计算，为14倍左右。各核心段获得加速效果不同，一方面是由于计算特征的差异导致与异构计算架构的可适应性不同，如计算密集型的程序相对于访存和通讯操作较多的程序更容易获得较好的加速效果。另一方面，优化效果不仅依赖于原有代码结构，同时与采用的优化思路和技术方法也密切相关。

将各优化后的核心段进行集成，形成BCCAGCM模式众核集成版本，分别采用100核组(6500核)、150核组(9750核)、200核组(13000核)、250核组(16250核)、300核组(19500核)和400核组(26000核)并行规模进行测试，由图 8可知，在各并行规模，采用从核对模式整体计算进行加速效果比较稳定，基本保持在1.9倍左右。

图 9为算例在不同并行规模加速比曲线，包括模式整体采用主核进行计算、模式整体采用主从核异构计算、模式物理过程采用主从核异构计算和模式动力框架采用主从核异构计算的加速比曲线。由图 9可知，模式采用主核、主从核进行计算时的加速趋势基本一致，在9750, 13000核可获得一定加速，在13000核时，加速比达到1.7左右。在16250, 19500核没有明显加速效果，这主要是因为模式在200~300 MPI进程(1个核组对应1个MPI进程)计算负载不平衡导致。26000核并行规模加速效果较为明显，加速比达到2以上。其中，采用异构计算的动力和物理部分相比较而言，现有优化版本，动力计算加速比达到2.8，效果优于物理计算。

与加速比对应，图 10给出算例在不同并行规模的并行效率，模式采用主核、主从核进行计算时并行效率变化趋势基本一致。在9750核规模，物理计算效率达到90%左右，其他均在87%左右。在13000核，动力计算略有上升，其他均呈下降趋势。在16250核、19500核规模，计算效率均呈下降趋势。在26000核规模，均有上升，动力计算并行效率约70%，其他约为57%。

由于气候模拟需要月、季节、年、十年甚至百年时间尺度的长期积分，对气候模式的计算精度有很高的要求。为了分析验证众核加速对模式模拟结果精度的影响，分别采用主核、主从核异构计算方式，设计两个平行数值模拟试验，两个试验采用相同的外强迫场、初值场，初始积分时间均从1978年1月1日开始，分别连续积分2年，输出月平均的物理量场。优化试验(主从核异构计算)和非优化试验(主核计算)模拟结果相比较, 两个试验都能合理再现全球年平均地表气温的空间分布特征，并且两个试验的全场最大值(305.71 K)、最小值(220.07 K)以及全球平均(287.49 K)完全一致。从差值场看，优化试验相对非优化试验的平均偏差为0，均方根误差也为0，两个试验结果一致。另外，两个试验模拟的冬季和夏季平均的地表气温气候场结果也一致(图略)。对比优化试验和非优化试验模拟的长期平均的年平均、冬季和夏季平均的34个变量的全球平均结果发现，两个试验模拟这些变量结果差值和均方根误差值均为0，即二者一致^[28-30]。

一般情况下，由于申威芯片主从核结构和数学库的差异，采用从核加速的计算结果与主核计算结果相比，会略有不同。但考虑到气候模式需要进行长时间积分，保持计算结果的精确性很重要，因此，关闭了一些优化功能，慎重选择调用的数学库，在损失一定计算效率的情况下，避免了计算误差，保证了此次研究工作计算结果的一致性。

4 小结与讨论

本文基于神威·太湖之光国产异构众核高性能计算系统对BCCAGCM模式进行了移植与优化，通过对模式代码进行性能分析，确定计算核心段，采用OpenACC方法对核心段代码进行加速，并完成众核版本集成，对长期积分结果进行合理性检验。研究表明：

1) 将BCCAGCM模式计算耗时占比较高的几个核心段，通过优化访存方式、优化数据传输策略、调整循环结构等方法，充分利用LDM空间，增加并行度。优化后，各核心段计算效率基本达到未优化计算的3倍左右，最高约14倍。

2) 对各优化后的模式核心段代码进行集成测试，模式整体计算加速是未优化计算的1.9倍左右，最大并行规模达到26000核，并行效率为57%。对2年模拟结果进行分析，众核集成版本可正确、稳定运行，计算误差合理。

3) 实践证明，BCCAGCM模式具有代码量大、计算逻辑复杂的特征，通过对代码并行算法进行深入剖析，结合异构硬件体系架构，对代码进行重构改写，可获得一定加速效果，为模式代码在异构计算平台的移植积累了技术经验。首次实现了我国自主研发的气象数值模式在我国自主研发的高性能计算系统大规模并行稳定运行。

由于本文采用了与模式开发不同的计算平台，优化过程中，检测出模式原有代码部分隐含问题，对其进行修订和改进，提高了模式整体健壮性和可移植性。对模式动力框架和物理过程的更细粒度并行优化方法和思路部分可应用于现有业务模式的计算优化，提升其计算效率和可扩展性，对业务发展有一定的贡献。同时，考虑程序在异构计算系统的性能优化是一个持续渐进、循环迭代的过程，可获得的性能提升很大程度上与投入研发的力量密切相关。对于已经相对稳定成熟的业务模式，通过持续的技术投入和积累，更深入的程序优化改造，实现在国产众核平台优化运行是很有希望的。

本文将BCCAGCM全球大气环流模式移植到全国产异构众核计算系统，对数值预报预测研究由模式到计算系统的全国产化应用进行了有益的探索，是气象数值模式对新型加速体系架构的全新尝试，取得一定成效，为BCCAGCM大规模并行复杂科学计算程序异构计算提供了一个基础版本。但是由于目前只对部分核心段进行优化，模式整体加速效果还有待进一步提高，未来将分别从动力框架和物理过程两方面顶层算法设计着手，对代码重构策略进一步优化，促进算法和硬件体系架构的适应调整，并对模式在异构系统的可扩展性开展深入研究和调优。另一方面，充分考虑高性能计算系统混合架构与编程模型对数值模式软件设计的影响，基于NVIDIA GPU，Intel Xeon Phi等新型架构，开展模式的移植优化工作，为未来高分辨率数值模拟工作的开展提供技术支撑。